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ABSTRACT 


SRRIT  is  a FORTRAN  program  to  calculate  an  approximate  orthonormal  basis 
for  a dominant  invariant  subspace  of  a real  matrix  A.  Specifically, 


given  an  inte 


$ger  m, 


SRRIT  attempts  to  compute  a matrix  Q with  m 


orthonormal  columns  and  real  quasi-triangular  matrix  T of  order  m 
such  that  the  equation 

AQ  = QT 

is  satisfied  up  to  a tolerence  specified  by  the  user.  The  eigenvalues 
of  T are  approximations  to  the  m largest  eigenvalues  of  A,  and  the 
columns  of  Q span  the  invariant  subspace  corresponding  to  those  eigen- 
values. SRRIT  references  A only  through  a user  provided  subroutine  to 
form  the  product  AQ;  hence  it  is  suitable  for  large  sparse  problems. 

^ 

*This  work  was  supported  in  part  by  the  Office  of  Naval  Research  under 
Contract  No.  N 00014-76-C-0391. 
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SRRIT  - A FORTRAN  Subroutine 
to  Calculate  the  Dominant  Invariant  Subspaces 
of  a Real  Matrix 

G.  W.  Stewart 


DESCRIPTION 


1.  Introduction 

The  program  described  in  this  paper  is  designed  primarily  to  solve 
eigenvalue  problems  involving  large,  sparse,  real  matrices.  The  programs 
attempt  to  calculate  a set  of  the  largest  eigenvalues  of  the  matrix  in 
question.  In  addition  they  calculate  a canonical  orthonormal  basis  for 


the  invariant  subspace  spanned  by  the  eigenvectors  and  principal  vectors 
corresponding  to  the  set  of  eigenvalues.  No  explicit  representation  of 
the  matrix  is  required;  instead  the  user  furnishes  a subroutine  to  cal- 
culate the  product  of  the  matrix  with  a vector. 

Since  the  programs  do  not  produce  a set  of  eigenvectors  corresponding 
to  the  eigenvalues  computed,  it  is  appropriate  to  begin  with  a mathematical 
description  of  what  is  actually  computed  and  how  the  user  may  obtain  eigen- 
vectors from  this  output  if  he  so  desires.  Let  A be  a matrix  of  order 
n with  eigenvalues  X..,\7 X ordered  so  that 


ACCESSION  for 


An  invariant  subspace  of  A is  any  subspace  Q.  for  which 
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i.e.,  the  subspace  is  transformed  into  itself  by  the  matrix  A. 

If  Q.  is  an  invariant  subspace  of  A and  the  columns  of 
Q = (q1,q2,...,qm)  form  a basis  for  Qj  then  Aq^  € Q,  and  hence  Aqi 
can  be  expressed  as  a linear  combination  of  the  columns  of  Q;  i.e., 
there  is  an  m-vector  t^  such  that  AQ  = Qt^.  Setting 

we  have  the  relation 

(1.1)  AQ  = QT  . 

In  fact  the  matrix  T is  just  the  representation  of  the  matrix  A in 
the  subspace  Q with  respect  to  the  basis  Q. 

If  x is  an  eigenvector  of  T corresponding  to  the  eigenvalue  X, 
then  it  follows  from  (1.1)  and  the  relation  Tx  = Xx  that 

(1.2)  A(Qx)  = \(Qx)  , 

so  that  Qx  is  an  eigenvector  of  A corresponding  to  the  eigenvalue  X. 

Thus  the  eigenvalues  of  T are  also  eigenvalues  of  A.  Conversely  if 

X , X.  ,...,X.  are  any  m eigenvalues  of  A that  are  distinct  from  the 
X1  *2 

other  n-m  eigenvalues,  then  there  is  a unique  invariant  subspace  of  dimen- 
sion m corresponding  to  these  eigenvalues;  i.e.,  the  eigenvalues  of  T 

in  (1.1)  are  precisely  X.  ,X.  ,...,X.  . 

h x2  "Sn 

If  |X.|  > |X.+1|,  then  there  is  a unique  dominant  invariant  subspace 
corresponding  to  X1,X2 Xi>  When  and  <^+1  exist,  c 
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The  subroutine  SRRIT  attempts  to  compute  a nested  sequence  of  orthonormal 
bases  of  . . . , (^.  Specifically,  if  all  goes  well,  the  subroutine  pro- 

duces a matrix  Q with  orthonormal  columns  having  the  property  that  if 
'V  > IW  then  q^,q2* . . . ,q^  span  (2^. 

The  case  where  and  are  a complex  conjugate  pair,  and 

hence  |X^_jJ  = |\jj»  is  treated  as  follows.  The  matrix  Q is  calculated 
so  that  the  matrix  T in  (1.1)  is  quasi-triangular;  i.e.,  T is  block 
triangular  with  lxl  and  2x2  blocks  on  its  diagonal.  The  structure  of 
a typical  quasi-triangular  matrix  is  illustrated  below  for  m = 6: 
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The  lxl  blocks  of  T contain  the  real  eigenvalues  of  A and  the  2x2 
blocks  contain  conjugate  pairs  of  complex  eigenvalues.  This  arrangement 
enables  us  to  work  entirely  with  real  numbers,  even  when  seme  of  the  eigen- 
values of  T are  complex.  The  existence  of  such  a decomposition  is  a 
consequence  of  Schur's  theorem  (see  [ 9 ]). 

The  eigenvalues  of  the  matrix  T computed  by  the  program  appear  in 
descending  order  of  magnitude  along  its  diagonal.  For  fixed  i let 
Q*1  = (q^»cl2,”*’qi^  311(1  let  T'*  be  the  leading  principal  submatrix  of  T 
of  order  i.  Then  if  the  i-th  diagonal  entry  of  T does  not  begin  a 2 x 2 
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block,  we  have 

AQ!i  = Q^T  1 . 

Thus  the  first  i columns  of  Q span  the  invariant  subspace  correspond- 
ing to  the  first  i eigenvalues  of  T.  When  |X^|  > |X^  | this  is  the 
unique  dominant  invariant  subspace  CX  • When  |X^|  = |X^+^ | the  columns 
of  span  a dominant  invariant  subspace;  but  is  is  not  unique,  since 

there  is  no  telling  which  comes  first,  X^  or  X.+^. 

Any  manipulations  of  A within  the  subspace  £ corresponding  to  Q 
can  be  accomplished  by  manipulating  the  matrix  T.  For  example, 

A*Q  - QTk  , 

so  that  if  f(A)  is  any  function  defined  by  a power  series,  we  have 

f (A)Q  - Qf (T)  . 

If  the  spectrum  of  A that  is  not  associated  with  Q is  negligible, 
considerable  work  can  be  saved  by  working  with  the  genera1 ly  much  smaller 
matrix  T in  the  coordinate  system  defined  by  Q.  If  explicit  eigen- 
vectors are  desired,  they  may  be  obtained  by  evaluating  the  eigenvectors 
of  T and  applying  (1.2).  The  programs  hqr2  in  [12]  and  HQR1  in  [7] 
will  evaluate  the  eigenvectors  of  a quasi-triangular  matrix. 


! 

I 
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2.  Usage 

SRRIT  is  a FORTRAN  subroutine  to  calculate  the  basis  for  described 
in  Section  1.  The  calling  sequence  for  SFFIT  is 

CALL  SRRIT (Q,AQ,ATQ,N,NV,M,EPS,MAXIT, START, T,ER, El, 

TYPE,RSD,RSDX) 

with  (.starred  parameters  are  altered  by  the  subroutine) 

*Q(N,M)  A real  array  that  on  return  contains  the  approximation 
to  Q.  Initially  Q may  contain  a starting  approxima- 
tion (cf.  START). 

*AQ(N,M)  A real  array  that  on  return  contains  the  product  AQ. 

ATQ  The  name  of  a FORTRAN  subroutine  that  computes  the  product 

AQ.  For  details  see  below. 

N The  order  of  A. 

*NV  The  number  of  vectors  to  compute.  On  return,  NV  contains 

the  number  of  columns  of  Q that  have  converged. 

M The  number  of  columns  of  Q.  M must  be  greater  than  or 

equal  to  NV. 

EPS  A convergence  criterion. 

MAXIT  An  integer  containing  the  maximum  number  of  iterations  to 

perform. 

START  An  integer  that  tells  the  initial  status  of  Q.  If 

START  < 0,  a starting  approximation  is  to  be  generated 
randomly.  If  START  > 0,  then  Q initially  contains  a 
starting  approximation;  and  if  START  > 1,  then  the  columns 
of  Q are  assumed  to  be  orthonormal. 


*T(M,M)  A real  array  that  on  return  contains  the  approximation  to 
the  matrix  T of  (1 . . 


r 

i , 
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*TYPE(M)  An  integer  array  whose  i-th  entry  on  return  is 

0 if  the  i-th  eigenvalue  is  real 

1 if  the  i-th  eigenvalue  is  the  first  of  a con- 
jugate pair  of  complex  eigenvalues 

2 if  the  i-th  eigenvalue  is  the  second  of  a conju- 
gate pair  of  complex  eigenvalues 

-1  if  the  i-th  eigenvalue  was  not  correctly 
determined 

*RSD(M)  A real  array  whose  i-th  entry  is  the  2-norm  of  the 

residual  associated  with  the  i-th  column  of  Q [cf.  (3.2) 
below] . 

*RSDX(M)  An  integer  array  whose  i-th  entry  is  the  iteration  at 
which  the  i-th  entry  of  RSD  was  computed. 

The  dimensions  in  the  parameter  descriptions  are  the  smallest  for 
which  the  program  will  work.  In  the  program  listed  here  they  are  set  for 
values  of  N up  to  five  hundred  and  M up  to  ten.  To  accomodate  larger 
problems,  change  the  dimension  500  to  the  largest  expected  value  of  N and 
the  dimension  10  to  the  largest  expected  value  of  M throughout  SRRIT  and 
its  auxiliary  subroutines  (n.b.  this  includes  the  dimension  information 
in  the  subroutine  calls  in  SRRSTP). 

The  user  may  furnish  a starting  approximation  to  the  matrix  Q in 
the  array  Q.  Actually  all  that  is  required  is  a set  of  vectors  whose 
column  space  approximates  If  such  a starting  approximation  is  fur- 

nished, the  parameter  START  should  be  set  greater  than  or  equal  to  zero. 

If  the  starting  vectors  are  orthonormal,  the  parameter  START  should  be  set 
positive.  If  START  is  negative,  Q is  initialized  with  random  numbers 
and  orthogonalized  to  provide  the  starting  approximation. 

The  user  is  required  to  furnish  a subroutine  to  calculate  the  product 


i 


AQ.  The  calling  sequence  for  this  subroutine  is 


with 
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CALL  ATQ(Q,AQ,L1,L2) 


Q(N,M) 


A real  array  containing  the  matrix  Q. 


AQ(N,M)  A real  array.  On  return  columns  LI  through  L2  of 
AQ  should  contain  the  product  of  the  matrix  A with 
columns  LI  through  L2  of  Q. 

LI  ) Integers  which  specify  which  columns  of  Q to  multiply 

L2J  by  the  matrix  A. 

A call  to  ATQ  causes  the  iteration  counter  to  be  increased  by  one,  so  that 
the  parameter  N1AXIT  is  effectively  a limit  on  the  number  of  calls  to  ATQ. 

The  convergence  criterion  is  described  in  detail  in  Sections  3 and  4. 
Essentially  the  matrices  Q and  T calculated  by  the  program  will  satisfy 


(2.1) 


(A.E)QINV  = qinvtiw 


where  NV  (on  return)  is  the  number  of  columns  that  have  converged  and 
E is  of  order  EPS.  From  this  it  is  seen  that  EPS  should  be  small 
compared  with  A.  The  criterion  insures  that  the  well-conditioned  eigen- 
values of  A will  be  calculated  accurately,  and  the  well -conditioned 
eigenvectors  can  be  calculated  accurately  from  Q and  T. 

The  rate  of  convergence  cf  the  i-th  column  of  Q depends  on  the  ratio 
|Xj^j/X.|.  For  this  reason  it  may  be  desirable  to  take  the  number  of  columns 
M of  Q to  be  greater  than  the  number  of  columns  NV  that  one  desires  to' 
compute.  For  example,  if  the  eigenvalues  of  A are  1.0,  0.9,  0.5,...  it 
will  pay  to  take  M = 2,  even  if  only  the  eigenvector  corresponding  to  1.0 
is  desired. 
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Since  SRRIT  is  designed  primarily  to  calculate  the  largest  eigenvalues 
of  a large  matrix,  no  provisions  have  been  made  to  handle  zero  eigenvalues. 
In  particular,  zero  eigenvalues  can  cause  the  program  to  stop  in  the  auxi- 
liary subroutine  0RTH. 

SRRIT  is  supported  by  a number  of  auxiliary  subroutines  (SRRSTP, 
RESID,GR0UP,0RTH,C$ND,RAND0M)  which  are  described  in  Section  5.  It  also 
requires  the  EISPACK  subroutines  0RTHES  and  0RTRAN  [7],  and  the  subroutines 
HQR3,  EXCHNG,  SPLIT,  and  QRSTEP  [11] . 

SRRIT  can  be  used  as  a black  box.  As  such  the  first  NV  vectors  it  re- 
turns will  satisfy  (2.1),  although  not  as  many  vectors  as  the  user  requests 
need  have  converged  by  the  time  MAXIT  is  reached.  However,  the  construction 
of  the  program  has  involved  a number  of  arbitrary  decisions.  Although  the 
author  has  attempted  to  make  such  decisions  in  a reasonable  manner,  it 
is  too  much  to  expect  that  the  program  will  perform  efficiently  on  all 
distributions  of  eigenvalues.  Consequently  the  program  has  been  written 
in  such  a way  that  it  can  be  easily  modified  by  someone  who  is  familiar 
with  its  details.  The  purpose  of  the  next  three  sections  is  to  provide 
the  interested  user  with  these  details. 


3.  Method 
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The  Schur  vectors  Q of  A are  computed  by  a variant  of  simultaneous 
iteration,  which  is  a generalization  of  the  power  method  for  finding  the 
dominant  eigenvector  of  a matrix.  The  method  has  an  extensive  literature 
[1,2, 3, 5, 8],  and  Rutishauser  [6]  has  published  a program  for  symmetric 
matrices  from  which  many  of  the  devices  in  SRRIT  have  been  drawn.  The 
method  about  to  be  described  has  been  analyzed  in  [10]. 

The  iteration  for  computing  Q may  be  described  briefly  as  follows. 
Start  with  an  n x m matrix  Qq  having  orthonormal  columns.  Given  Q , 
form  Qv+1  according  to  the  formula 

Vi ' - 

where  **1  is  either  an  identity  matrix  or  an  upper  triangular  matrix 
chosen  to  make  the  columns  of  C>v+1  orthonormal  (just  how  often  such  an 
orthogonalization  should  be  performed  will  be  discussed  below).  If 
l\J  > I \n+l I ’ then  under  mi*d  restrictions  on  Qq  the  column  space  of 
Qv  approaches  C^. 

The  individual  columns  of  Qv  will  in  general  approach  the  correspond- 
ing columns  of  the  matrix  Q defined  in  Section  1;  however  the  rate  of 
convergence  of  the  i-th  column  is  proportional  to  max  1|v, 

|^+jAjJV}  and  may  be  intolerably  slow.  The  process  may  be  accelerated  by 
the  occasional  application  of  a "Schur-Rayleigh-Ritz  step"  (from  which 
SRRIT  derives  its  name),  which  will  now  be  described.  Start  with  Q just 
after  an  orthogonalization  step,  so  that  = I.  Form  the  matrix 
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B = QTAQ  , 

V V 

and  reduce  it  to  ordered  quasi- triangular  form  by  an  orthogonal 
similarity  transformation  Y^: 

(3.1)  YTB  Y = T . 

V V V V 

Finally  overvvTite  Q with  Q Y . 

The  matrices  formed  in  this  way  have  the  following  property. 

If  > |X^|  > |\i+1|,  then  under  mild  restrictions  on  Qq  the  i-th 

column  of  Qv  will  converge  to  the  i-th  column  q^  of  Q at  a 

rate  proportional  to  |Xm+1Ai|v.  Thus  not  only  is  the  convergence  accelerated, 
but  the  first  columns  of  Qv  tend  to  converge  faster  than  the  later  ones. 

A number  of  practical  questions  remain  to  be  answered. 

1.  How  should  one  determine  when  a column  of  Qv  has  converged? 

2.  Can  one  take  advantage  of  the  early  convergence  of  some  of  the 
columns  of  Qv  to  save  computations? 

3.  How  often  should  one  orthogonalize  the  columns  of  the  Q^? 

4.  How  often  should  one  perform  the  SRR  acceleration  described  above? 
Here  we  shall  merely  outline  the  answers  to  these  questions.  The  details 
will  be  given  in  the  discussion  of  SRRIT. 

1.  Convergence.  If  |X^j|  = IVj  or  |X^|  ■ |Xj+^|,  the  i-th  column 
of  Q is  not  uniquely  determined;  and  when  |X^|  is  close  to  |X^ 
or  |X^|,  the  i-th  column  cannot  be  computed  accurately.  Thus  a convergence 
criterion  based  on  the  i-th  column  q^  of  Qv  becoming  stationary  is 
likely  to  fail  when  A has  equimodular  eigenvalues.  Accordingly  we  have 
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adopted  a different  criterion  which  amounts  to  requiring  that  the  relation 

(1.1)  is  almost  satisfied.  Specifically,  let  t : denote  the  i-th 
column  of  in  (3.1).  Then  the  i-th  column  of  the  Qv  produced  by 
the  SRR  step  is  said  to  have  converged  if  the  2-norm  (see  [9]  for  a defi- 
nition) of  the  residual  vector 

(3.2)  r M - AqM  ..  QvtW 

is  less  than  some  prescribed  tolerance. 

If  this  criterion  is  satisfied  for  each  column  of  Q , then  the  resi- 

v 

dual  matrix 

R = AQ  - Q T 

V V V 

will  be  small.  This  in  turn  implies  that  there  is  a small  matrix 
E = -R  0T  such  that 

v VXV 

(A+F.  )Q  = Q T , 

so  that  Qv  and  T^  are  the  matrices  associated  with  the  slightly  per- 
turbed matrix  A + E^,  provided  only  that  some  small  eigenvalue  of  A + E^ 
has  not  by  happenstance  been  included  in  T^.  To  avoid  this  possibility  we 
group  nearly  equimodular  eigenvalues  together  and  require  that  their  average 
value  has  settled  down  before  testing  their  residuals.  In  addition  a group 
of  columns  is  tested  only  if  the  preceding  columns  have  all  converged. 

2.  Deflation.  The  theory  of  the  iteration  indicates  that  the  initial 
columns  of  the  Qv  will  converge  before  the  later  ones.  When  this  happens 
considerable  computation  can  be  saved  by  freezing  these  columns.  This  saves 


multiplying  the  frozen  columns  by  A,  orthogonal izing  them  when  + I, 

and  work  in  the  SRR  step. 

3.  Orthogonalization.  The  orthogonalization  of  the  columns  of  AQ^ 
is  a moderately  expensive  procedure  which  is  to  be  put  off  as  long  as 
possible.  The  danger  in  postponing  orthogonalization  is  that  cancellation 
of  significant  figures  can  occur  when  AQ^  is  finally  orthogonalized, 

as  it  must  be  just  before  an  SRR  step.  In  [10]  it  is  shown  that  one  can 
expect  no  more  than 

(3.3)  t * k log10  k(T) 

decimal  digits  to  cancel  after  k iterations  without  orthogonalization 
(here  k(T)  = ||T||||T  ^||  is  the  condition  number  of  T with  respect  to  inver- 
sion). The  relation  (3.3)  can  be  used  to  determine  the  number  of  iterations 
between  orthogonal izat ions . 

4.  SRR- steps.  The  SRR- step  described  above  does  not  actually  acceler- 
ate the  convergence  of  the  Qv;  rather  it  unscrambles  approximations  to  the 
columns  of  tliat  are  already  present  in  the  colimn  space  of  and 
orders  them  properly.  Therefore,  the  only  time  an  SRR  step  needs  to  be 
performed  is  when  it  is  expected  that  a column  has  converged.  Since  it  is 
known  from  the  theory  of  the  iteration  that  the  residuals  in  (3.2)  tend 
linearly  to  zero,  the  iteration  at  which  they  will  satisfy  the  convergence 
criterion  can  be  predicted  from  their  values  at  two  iterations.  As  with 
convergence,  this  prediction  is  done  in  groups  corresponding  to  nearly 
equLmodular  eigenvalues. 


4.  Details  of  SRRIT 


In  designing  SRRIT,  we  have  tried  to  make  it  easily  modifiable.  This 
has  been  done  in  two  ways.  First,  we  have  defined  a number  of  important 
control  parameters  and  given  them  values  at  the  beginning  of  the  program. 

The  knowledgeable  user  may  alter  these  values  to  improve  the  efficiency 
of  the  program  in  solving  particular  problems.  Second,  a number  of  important 
tasks  have  been  isolated  in  independent  subroutines.  This  should  make  it 
easy  to  modify  the  actual  structure  of  SRRIT,  should  the  user  decide  that 
such  radical  measures  are  necessary.  In  this  section  we  shall  describe 
SRRIT  in  some  detail,  specifying  the  action  of  the  control  parameters.  In 
the  next  section  we  shall  describe  the  supporting  subroutines. 

Here  follows  a list  of  the  control  parameters  with  their  initial  values 
and  a brief  description  of  their  functions. 


I NIT  (5) 

STPFAC  (2.0) 

ALPHA  (1.0) 
BETA  (1.1) 

GRPTOL  (0.001) 

CNVTOL  (0.001) 

ORTTOL  (2.0) 

SEED  (69) 


a number  of  initial  iterations  to  be  performed  at 
the  outset. 

a constant  used  to  determine  the  maximum  number  of 
iterations  before  the  next  SRR  step. 

parameters  used  in  predicting  when  the  next  residual 
will  converge. 

a tolerance  for  grouping  equimodular  eigenvalues 

a convergence  criterion  for  the  average  value  of  a 
cluster  of  equimodular  eigenvalues. 

the  nunber  of  decimal  digits  whose  loss  can  be 
tolerated  in  orthogonal ization  steps. 

a seed  for  the  random  nunber  generator  that  initializes  Q. 
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We  now  give  an  informal  description  of  SRRIT  as  it  appears  in  the 
ALGORITHM  section.  The  variable  L points  to  the  first  column  of  Q 
that  has  not  converged.  The  variable  IT  is  the  iteration  counter.  The 
variable  NXTSRR  is  the  iteration  at  which  the  next  SRR  step  is  to  take 
place,  and  the  variable  DORT  is  the  interval  between  orthogonalizations . 


srrit: 

1.  initialize  control  parameters 

2.  initialize 

1.  IT  = 0 

2.  L = 1 

3.  Initialize  Q as  prescribed  by  START 

3.  srri  loop 

1.  perform  an  SRR  step 

2.  confute  the  residuals 

3.  check  convergence,  resetting  L if  necessary 
4..  if  L > NV  or  IT  > MAXIT  then  leave  srr 

5.  calculate  KXTSRR 

6.  calculate  D0RT  and  NXT0RT 

7.  Q = AQ;  IT  = IT+1 

8.  orth:  loop  until  IT  = NXTSRR 

1.  power:  loop  until  IT  = NXT0RT 

1.  AQ  =T*Q 

2.  Q = AQ 

3.  IT  = IT+1 
end  power 

2.  orthogonalize  Q 

3.  NXT0RT  = min  (NXTSRR, IT+D0RT) 
end  orth 

end  srr 

4.  NV  * L-l 
end  srrit 

The  details  of  this  outline  are  as  follows  (the  numbers  correspond  to  the 
statements  in  the  algorithm). 

2.3.  If  START  < 0,  then  Q is  initialized  using  the  function  RANDOM. 
If  START  £ 0,  the  columns  of  Q are  orthogonal i zed  by  the  subroutine  0RTH. 


3.  This  is  the  main  loop  of  the  program.  Each  time  it  is  executed 
an  SRR  step  is  performed  and  convergence  is  tested. 

3.1.  The  SRR  step  is  performed  by  the  subroutine  SRRSTP,  which 
returns  the  new  Q and  A*Q,  as  well  as  T and  its  eigenvalues. 

3. 2.  The  residuals  are  computed  by  the  subroutine  RESID. 

3.3.  The  algorithm  for  determining  convergence  is  the  following. 

Starting  with  the  L-th  eigenvalue,  the  subroutine  GR0UP  is  called  to  deter- 
mine a group  of  nearly  equimodular  eigenvalues,  as  defined  by  the  parameter 
GRPT0L.  The  same  is  done  for  the  old  eigenvalues  from  the  last  SRR  step. 

If  the  groups  have  the  same  number  of  eigenvalues  and  the  average  value  of 
the  eigenvalues  has  settled  down  (CNVT0L) , then  the  residuals  are  averaged 
and  tested  against  EPS.  If  the  test  is  successful,  L is  increased  by 

the  number  in  the  group,  and  the  tests  are  repeated.  Otherwise  control  is 

passed  to  statement 

3.4.  where  the  two  termination  conditions  for  SRRIT  are  tested. 

3.5.  The  iteration  at  which  the  next  SRR-step  is  to  take  place 
(NXTSRR)  is  determined  as  follows.  NXTSRR  is  tentatively  set  equal  to 
STPFAC*IT.  If  the  number  of  eigenvalues  in  the  new  and  old  groups  corre- 
sponding to  the  next  set  of  unconverged  eigenvalues  is  the  same,  the  RMS 
average  of  the  norms  of  the  residuals  of  each  group  is  calculated 
(ARSD,  A0RSD) . If  ARSD  < EPS,  then  NXTSRR  = IT+1.  If  ARDS  > A0RSD,  then 
NXTSRR  = STPFAC*IT.  Otherwise 


NXTSRR  - min  (IT+ALPHA+BETA*DSRR, STPFAC*IT) 


v/here 


r 
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DSRR  - (0RSDX-RSDX)  . 

Finally  NXTSRR  is  constrained  to  be  less  than  or  equal  to  MAXIT. 

3.6.  The  interval  D0RT  between  orthogonal izat ions  is  computed  from 
(3.3): 


D0RT  = max  (l,0RIT0L/loglo  k(T)), 

where  the  condition  number  k(T)  is  calculated  by  the  function  C0ND. 

The  next  orthogonalization  occurs  at 

NXT0RT  = min  (IT+D0RT, NXTSRR)  . 

3.7.  Since  the  SRR  step  computes  a product  AQ,  the  iteration  count 
must  be  increased  and  AQ  placed  back  in  Q. 

3.8.  Loop  on  orthogonal izations. 

3.8.1.  Loop  overwriting  Q with  the  product  A*Q. 

4.  Set  NV  to  the  number  of  vectors  that  have  actually  converged  and 


return. 


5.  Auxiliary  Subroutines 

In  this  section  we  shall  describe  the  subroutines  called  by  SRRIT. 

Some  of  these  subroutines  have  been  coded  in  greater  generality  than  is 
strictly  required  by  SRRIT  in  order  to  make  the  program  easily  modifiable 
by  the  user. 

SRRSTP(Q,AQ,ATQ,L,M,N,T,ER,EI,TYPE,0ER,0EI,0TYPE)  . 

This  subroutine  performs  an  SRR  step  on  columns  L through  M of 
Q.  After  forming  AQ  and  T = QT (AQ) , the  routine  calls  0RTHES,  0RTRAN, 
and  HQR3  to  reduce  T to  ordered  quasi-triangular  form.  The  triangularizing 
transformation  is  postmultiplied  into  Q and  AQ.  The  eigenvalues  from  the 
last  step  are  stored  in  the  arrays  0ER,  0EI  and  0IYPE,  and  the  new  eigenvalues 

are  placed  in  the  arrays  ER,  El , and  TYPE. 

RESID(Q,AQ,T,RSD,RSDX,0RSD,0RSDX,L1,L2,M,N,IT,TYPE)  . 

This  subroutine  computes  the  norm  of  the  residuals  (3.2)  for  columns 
LI  through  L2  of  Q.  The  old  residuals  and  their  iteration  numbers  are 
saved  in  the  arrays  0RSD  and  0RSDX.  The  I-th  entry  of  the  array  RSDX  is 
set  to  IT  depending  on  whether  or  not  TYPE (I)  > 0.  For  a complex  pair  of 
eigenvalues,  the  RMS  average  of  the  norms  of  their  two  residuals  is  re- 
turned. 

GR0UP (ER, El ,TYPE , GRPT0L , L,M,N,NGRP , CTR, AE) 

This  subroutine  locates  a group  of  approximately  equimodular  eigenvalues 
^L’^L+l’ ' * ’ ,XL+NGRP-1  T^e  eigenvalues  so  grouped  satisfy 


^ -CTR 


< GRPTOL*CTR  (i-L,L+l, . . . ,L+NGRP-1)  . 


The  mean  of  the  group  is  returned  in  AE. 

0RTH(AQ,Q,L,M,N) 

For  J = L,L+1,...,M  this  subroutine  orthogonalizes  the  J-th  column 
of  AQ  with  respect  to  columns  1,2, . . . ,L-1  of  Q and  columns 
L,L+1, . . . ,J-1  of  AQ.  The  results  are  returned  in  Q.  The  method  used 
is  the  modified  Gram-Schmidt  method  with  reorgonalization.  No  more  than 
NTRY  reorthogonalizations  are  performed,  after  which  the  routine  executes 
a ST0P.  The  routine  will  also  stop  if  any  colunn  becomes  zero. 

RAND0MCSEED) 

This  function  subprogram  returns  a floating-point  pseudo-random  number 
between  0 and  1.  It  is  used  to  initialize  Q. 


6.  Numerical  Examples 
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The  program  described  above  has  been  tested  on  a number  of  problems. 
In  this  section  we  give  two  examples  that  illustrate  the  flexibility  of 
the  method  and  its  ability  to  deal  with  equimodular  eigenvalues. 

The  first  example  is  a random  walk  on  an  (n+1)  x (n+1)  triangular 
grid,  which  is  illustrated  below  for  n = 6. 

6 . 

5 . 

4 . . . 

3 . . . . 

2 

1 .....  . 

0 

V/h  0 1 2 3 4 5 6 

The  points  of  the  grid  are  labelled  (v,h)  (v«0,...,n-h;  h-0,...,n).  From 
the  point  (v,h) , a transition  may  take  place  to  one  of  the  forr  adjacent 
points  (v+l,h),  (v,h+l),  (v-l,h),  and  (v,h-l).  The  probability  of  jumping 
to  (v-l,h)  or  (v„h-l)  is 

(6*1)  pd(v,h)  ■ (v+h)/n 

with  the  probability  being  split  equally  between  the  two  points  when  both 
are  on  the  grid.  The  probability  of  junping  to  (v+l,h)  or  (v,h+l)  is 

(6*2)  pu(v,h)  - 1 - pd  (v,h) 

with  the  probability  again  being  split  when  both  points  are  on  the  grid. 


i 
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If  the  (n+l)(n+2)/2  nodes  (v,h)  are  numbered  1,2, . . . , (n+1) (n+2)/2 
in  some  fashion,  then  the  random  walk  can  be  expressed  as  a finite  Markov 
chain  whose  transition  matrix  A consists  of  the  probabilities  a^  of 
jumping  from  node  j to  node  i (A  is  actually  the  transpose  of  the 
usual  transition  matrix;  see  [4]).  To  calculate  the  i-th  element  of  the 
vector  Aq  one  need  only  regard  the  components  of  q as  the  average 
nunber  of  individuals  at  the  nodes  of  the  grid  and  use  the  probabilities 

(6.1)  and  (6.2)  to  calculate  how  many  individuals  will  be  at  node  i 
after  the  next  transition. 

We  are  interested  in  the  steady  state  probabilities  of  the  chain, 
which  is  ordinarily  the  appropriately  scaled  eigenvector  corresponding  to 
the  eigenvalue  unity.  However,  if  we  number  the  diagonals  on  the  grid 
that  are  parallel  to  the  hypotenuse  by  0,l,2,...,n,  then  an  individual 
on  an  even  diagonal  can  only  jimp  to  an  odd  diagonal,  and  vice  versa.  This 
means  that  the  chain  is  cyclic  with  period  two.  Computationally  it  means 
that  A has  an  eigenvalue  of  -1  as  well  as  +1. 

To  run  the  probl'an  on  SRRIT,  the  nodes  of  the  grid  were  matched  with 
the  components  of  the  vector  q in  the  order  (0,0) , (1,0) , . . . , (n,0) , (0,1) , 

(1.1) ,...,(n,l),(0,2) The  subroutine  that  computes  AQ  is  listed 

in  the  appendix.  Note  that  the  matrix  A is  never  explicitly  used;  all 
computations  are  done  in  terms  of  the  transition  probabilities  (6.1)  and 

(6.2) .  The  use  of  a common  block  to  transmit  information  from  the  program 
that  called  SRRIT  is  typical. 
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The  problem  was  run  for  a 30  * 30  grid  which  means  N = 496.  We 
took  M = 6,  NV  = 4,  and  EPS  = 10  The  results  for  each  iteration  in 
which  an  SRR  step  was  performed  are  summarized  in  the  following  table. 
The  variables  ER  and  El  are  the  real  and  imaginary  parts  of  the 
eigenvalues  and  RSD  is  the  norm  of  the  corresponding  residual.  CTR 
is  the  center  of  the  current  convergence  cluster,  AE  is  the  average 
value  of  the  eigenvalues  in  the  cluster,  and  ARSD  is  the  RMS  average 
of  the  residuals.  DSRR  is  the  number  of  iterations  to  the  next  SRR 
step  and  D0RT  is  the  number  to  the  next  orthogonalization. 
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IT  = 0 


ER 

.9457+00 

-.9096-01 

-.4841-01 

.3469-01 

.3469-01 

-.1921-01 

El 

0 

0 

0 

.1617-01 

-.1617-01 

0 

RSD 

CTR 

AE 

ARSD 

.38+00 

.9457+00 

.9457+00 

.38+00 

.60+00 

.61+00 

DSRR=5 

. 59+00 

D0RT=1 

.59+00 

.63+00 

IT  = 5 

ER  , 

.1012+01 

-.3912+00 

.2400+00 

-.1800+00 

.1371+00 

.3517-01 

El 

0 

0 

0 

0 

0 

0 

RSD 

.19+00 

.84+00 

. 93+00 

.89+00 

.91+00 

.93+00 

CTR 

.1012+01 

AE 

.1012+01 

DSRR=5 

D0RT=1 

ARSD 

.19+00 

IT  = 10 

ER 

.1017+01 

-.5987+00 

-.3499+00 

.3251+00 

.9255-01 

.5706-01 

El 

0 

0 

0 

0 

0 

0 

RSD 

.12+00 

.75+00 

.89+00 

.92+00 

. 95+00 

. 95+00 

CTR 

.1017+01 

AE 

.1017+01 

DSRR=10 

D0RT=1 

ARSD 

.12+00 

IT  = 20 

ER 

.1009+01 

-.8751+00 

.5175+00 

-.5124+00 

.3747+00 

-.1485+00 

El 

0 

0 

0 

0 

0 

0 

RSD 

. 58=01 

.46+00 

.82+00 

.84+00 

.88+00 

. 94+00 

CTR 

.1009+01 

AE 

.1009+01 

DSSR=20 

D0RT=2 

ARSD 

.58-01 

IT  = 40 


ER 

.1001+01 

-.9843+00 

.9195+00 

-.9144+00 

.7946+00 

-.5166+00 

El 

0 

0 

0 

0 

0 

0 

RSD 

CTR 

AE 

ARSD 

.23-01 

.1001+01 

*1001+01 

.23-01 

.14+00 

.37+00 

DSSR=40 

.40+00 

D0RT=1 

. 55+00 

.95+00 
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IT  = 80 


ER 

.1000+01  -.9998+00 

.9935+00  -.9934+00 

.8734+00 

El 

0 0 

0 0 

0 

RSD 

.74-02  .29-01 

.36-01  .78-02 

.43+00 

CTR 

.1000+01 

AE 

.1991-93 

DSRR=80  D0RT=5 

ARSD 

.21-01 

IT  = 160 

ER 

-.1000+01  .1000+01 

.9935+00  -.9935+00 

.9470+00 

El 

0 0 

0 0 

0 

RSD 

.56-03  .13-02 

.38-03  .70-03 

.23-03 

CTR 

.1000+01 

AE 

-.1304-04 

DSRR=135  D0RT=2 

ARSD 

.10-02 

IT  = 295 

ER 

-.1000+01  .1000+01 

.9935+00  -.9935+00 

.9755+00 

El 

0 0 

0 0 

0 

RSD 

.30-04  .37-05 

.13-05  .12-03 

.84-02 

CTR 

.1000+01 

.9935+00 

AE 

-.1863-06 

.1080-06 

DSRR=30 

ARSD 

.21-04 

.83-04 

IT  = 325 

ER 

-.1000+01  .1000+01 

.9935+00  -.9935+00 

.9755+00 

El 

0 0 

0 0 

0 

RSD 

.70-05  .82-06 

.35-06  .34-04 

.39-02 

CTR 

.1000+01 

.9935+00 

AE 

-.4470-07 

.7451-07 

DSRR=23 

ARSD 

.50-05 

.24-04 

IT  = 348 

ER 

.9935+00  -.9935+00 

.9755+00 

El 

0 0 

0 

RSD 

.12-06  .12-04 

.21-02 

CTR 

.9935+00 

AE 

.1118-06 

ARSD 

.88-05 

.2408+00 

0 

.92+00 


-.2138+00 

0 

.94-03 


-.9738+00 

0 

.57-01 

D0RT=3O 


-.9751+00 

0 

.26-01 


D0RT=23 


-.9754+00 

0 

.15-01 
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The  course  of  the  iteration  is  unexceptionable.  The  program  doubles 
the  interval  between  SRR  steps  until  it  can  predict  convergence  of  the 
first  cluster  corresponding  to  the  eigenvalues  ±1.  The  first  prediction 
falls  slightly  short,  but  the  second  gets  it.  After  a third  prediction 
the  program  terminates  on  the  convergence  of  the  second  group  of  two  eigen- 
values. 

It  should  be  noted  that  the  eigenvalue  -1  has  appeared  as  the 
dominant  one.  A transformation  bringing  the  eigenvector  corresponding  to 
1 can  be  obtained  by  calling  EXCHNG  of  [li ] to  interchange  the  eigenvalues 
1 and  -1  (however,  in  this  case  the  eigenvector  corresponding  to  1 is 
just  the  absolute  value  of  the  eigenvector  corresponding  to  -1). 

Without  actually  making  timing  runs,  it  is  difficult  to  predict  how 
much  work  is  entailed  in  finding  the  eigenvalues.  For  example,  runs  were 
made  with  M = 2,4, 6,8,  which  gave  the  following  table  of  iterations 
required  for  the  convergence  of  the  first  group  of  two  eigenvalues. 

m it  m*it 

2 1737  3474 

4 523  2092 

' . 6 325  1950 

8 188  1504 

As  predicted  by  the  convergence  theory,  the  number  of  iterations  decreases 
as  m increases.  However,  as  m increases  we  must  also  multiply  more 
columns  of  Q by  A,  and  for  this  particular  problem  the  number  m*it  is 
probably  a better  measure  of  the  amount  of  work  involved.  From  the  table  it 
is  seen  that  this  measure  is  also  decreasing,  although  less  dramatically  than 


~ 
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the  number  of  iterations.  This  of  course  does  not  include  the  overhead 
generated  by  SRRIT  itself,  which  increases  with  m and  may  be  considerable. 

The  second  example  shows  how  SRRIT  can  be  used  in  conjunction  with 
the  inverse  power  method  to  find  the  smallest  eigenvalues  of  a matrix. 
Consider  the  boundary  value  problem 

y"  + p2y  = 0 , 

(6.3)  y(0)  = 0 , 

y'(0)  + Yy'(l)  - 0,  0 < Y < 1 . 

The  eigenvalues  of  this  problem  are  easily  seen  to  be  given  by 

p.  * i cosh  ^ (-y  ^)  , 

which  are  complex.  The  following  table  lists  the  reciprocals  of  the  first 
eight  eigenvalues  for  y - 0.01. 


-2 

l -2i 

U 

lu  1 

0.01264  ± 0. 02313i 

. 02636 

0.004446  ± 0.007308 

.008544 

0.002895  ± 0.002204 

. 003638 

0.001740  ± 0.0008901 

.001954 

The  solution  of  (6.3)  can  be  approximated  by  finite  difference 
techniques  as  follows.  Let  y^  denote  the  approximate  solution  at  the 
point  x^  » i/(n+l)  (i*0,l,...,n+l).  Replacing  the  derivatives  in  (6.3) 

with  three  point  difference  operators,  we  obtain  the  following  generalized 
matrix  eigenvalue  problem  for  y - (y1,y2i..*»>rn+1)T: 


where 
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Ay  + M-^By  = 0 , 


A = 


-2  1 
1-2  1 
0 1-2 


0 


o 


0 

1 


1 0. . . 0 y -4y  3y 


2 

and  B = h diag  (1,1, . . . ,1,0) • We  may  recast  this  problem  in  the  form 


cy  = X y , 

M- 

where  C = A 1B. 

To  apply  SRRIT  to  this  problem,  we  must  be  able  to  compute  z = Cq  for 
any  vector  q.  This  can  be  done  by  solving  the  linear  system 


Az  = Bq  , 

which  is  easily  done  by  sparse  Gaussian  elimination. 

The  problem  was  run  for  n = 300  with  M = 6,  NV  = 4,  and  EPS  = 10’4. 
The  results  were  the  following: 


ER 

-.1525+00 

.1179-02 
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IT  = 0 

.1548-03 

.9887-04 

.9887-04 

.2577-04 

El 

0 

0 

0 

.5598-04 

-.5598-04 

0 

RSD 

.15+00 

.85-02 

.65-02 

.15-01 

.15-01 

.71-02 

CTR 

AE 

ARSD 

.1525+00 

-.1525+00 

.15+00 

DSRR=5 

D0RT=1 

IT  = 5 

ER 

-.1264-01 

-.1264-01 

.4438-02  .4438-02 

.3104-02 

.3104-02 

El 

.2313-01 

-.2313-01 

.7323-02  -.7323-02 

.2402-02 

-.2402-02 

RSD 

.85-07 

.85-07 

.81-05  .81-05 

.20-03 

.20-03 

CTR 

.2636-01 

AE 

.1264-01 

DSRR=5  D0RT=1 

ARSD 

85-07 

IT  = 10 

ER 

-.1264-01 

-.1244-01 

.4447-02  .4447-02 

.2909-02 

.2909-02 

El 

.2313-01 

-.2313-01 

.7308-02  -.7308-02 

.2204-02 

-.2204-02 

RSD 

.60-08 

.60-08 

.16-07  .15-07 

.93-05 

.93-05 

CTR 

.2636-01 

.8555-02 

AE 

-.1264-01 

.4447-02 

ARSD 

60-08 

.15-07 

Given  the  extremely  favorable  ratios  of  the  eigenvalues  in  Table  (6.4) 
--the  absolute  value  of  the  ratio  of  the  seventh  to  the  first  is  about  .075 
--it  is  not  surprising  that  the  iteration  converges  quickly.  Indeed  the 
only  thing  preventing  convergence  at  the  fifth  iteration  is  that  the  first 
eigenvalue  changed  from  real  in  the  first  iteration  to  complex  in  the  fifth. 
Thus  the  problem  is  hardly  a fair  test  of  machinery  of  SRRIT.  However,  it  is 
an  excellent  example  how  easy  it  is  to  apply  SRRIT  to  a problem  with  complex 
eigenvalues.  It  also  disposes  of  the  notion  that  large  eigenvalue  problems 
must  always  require  a large  amount  of  work  to  solve;  the  factor  that  limits 
the  size  of  n is  the  storage  available,  not  the  time  required  to  compute  Ax. 
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SUBROUTINE  SRRIT <Q,AQ ,ATQ ,N ,NV . M. EPS . MAXIT , START , T. 

1 ER>EI , TYPE, RSD,RSDX, WRITE) 

C 

C PARAMETERS  IN  THE  CALLING  SEQUENCE 
C 

INTEGER  N, NV,M,MAXIT, START, TYPE < 10 ) ,RSDX( 10) 

REAL  Q(500, 10) ,AQ<300, 10) ,EPS,T( 10, 10 > ,ER < 10 ) ,EI < 10 > . RSD< 10) 
LOGICAL  WRITE 
EXTERNAL  ATQ 
C 


C SRRIT  IS  A FORTRAN  SUBROUTINE  TO  COMPUTE  A NESTED  SEQUENCE 
C OF  ORTHONORMAL  BASES  FOR  THE  DOMINANT  INVARIANT  SUBSPACES  OF 

C A REAL  MATRIX  A OF  ORDER  N.  SPECIFICALLY,  THE  PROGRAM  RETURNS 

C AN  NXNV  MATRIX  Q WITH  ORTHONORMAL  COLUMNS  AND  AN  NV  X NV  MATRIX  T 

C SATISFYING 

C 

C A*Q  « Q*T  + O(EPS) . 

C 

C THE  MATRIX  T IS  QUASI -TRIANGULAR , THAT  IS  IT  IS  BLOCK 

C TRIANGULAR  UITH  1X1  AND  2X2  BLOCKS  ON  ITS  DIAGONAL.  THE 

C EIGENCALUES  IN  THE  1X1  BLOCKS  ARE  REAL.  THE 

C THE  EIGENVALUES  IN  THE  2X2  BLOCKS  ARE  COMPLEX  CONJUGATE 

C PAIRS.  THE  EIGENVALUES  E<1),  E(N>  ARE  ORDERED 

C SO  THAT 

C 

C ABS ( E ( 1 ) > .GE.  ABS (E < 2 > > . GE GE . ABS(E(NV>), 

C 

C AND  THESE  EIGENVALUES  APPROXIMATE  THE  LARGEST  EIGENVALUES 

C OF  A.  THESE  FACTS  HAVE  THE  FOLLOWING  CONSEQUENCES. 

C 

C 1.  IF  E(L>  .NE.  E(L+1 ) AND  E(L)  .NE.  CONJ ( E ( L+l ) ) , 

C THEN  COLUMNS  1»2»...,L  OF  Q FORM  AN  APPROXIMATE 

C BASIS  FOR  THE  INVARIANT  SUBSPACE  CORRESPONDING  TO 

C THE  L LARGEST  EIGENVALUES  OF  A.  THE  LXL  LEADING 

C PRINCIPAL  SUBMATRIX  OF  T IS  A REPRESENTATION  OF 

C A IN  THAT  SUBSPACE  WITH  RESPECT  TO  THE  BASIS  Q. 

C 

C 2.  IF  Z IS  AN  EIGENVECTOR  OF  T CORRESPONDING  TO  E, 

C THEN  Q*Z  IS  AN  APPROXIMATE  EIGENVECTOR  OF  A 

C CORRESPONDING  TO  E. 

C 

C THE  PROGRAM  ACTUALLY  ITERATES  WITH  AN  NXM  MATRIX  Q 
C AND  AN  MXM  MATRIX  T.  SINCE  THE  RATE  OF  CONVERGENCE 
C OF  THE  L-TH  COLUMN  OF  Q IS  ESSENTIALLY  LINEAR  WITH 

C RATIO  ABS <E(M+1)/E<L))»  IT  MAY  PAY  THE  USER  TO  SET  M 

C LARGER  THAN  THE  NUMBER,  NV,  OF  VECTORS  HE  WANTS  TO 

C COMPUTE . 

C 

C THE  USER  MUST  FURNISH  A SUBROUTINE  TO  COMPUTE  THE 
C PRODUCT  A*Q.  THE  CALLING  SEQUENCE  IS 

C 

C CALL  ATQ (Q, AO, LI ,L2) 

C 

C FOR  J-LlfLl+l, ...»L2  THE  PROGRAM  MUST  PLACE  THE  PRODUCT 
C A*Q(*,J>  IN  AQ(*,J>. 

C 

C THE  PARAMETERS  IN  THE  CALLING  SEQUENCE  OF  SRRIT  ARE 
C (STARRED  PARAMETERS  ARE  ALTERED  BY  THE  PROGRAM) 

C 

C *Q  AN  ARRAY  THAT  ON  RETURN  CONTAINS  THE 

C ORTHONORMAL  VECTORS  DESCRIBED  ABOVE.  INITIALLY 

C Q MAY  CONTAIN  A STARTING  APPROXIMATION 

C <CF . START  BELOW). 

C *AQ  AN  ARRAY  THAT  ON  RETURN  CONTAINS  THE  PRODUCT 
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MO. 

ATQ  THE  NAME  OF  A SUBROUTINE  TO  EVALUATE  THE 

PRODUCT  A*0. 

N THE  ORDER  OF  THE  HATRIX  A. 

«NV  THE  NUMBER  OF  VECTORS  TO  COMPUTE.  ON  RETURN > 

NV  CONTAINS  THE  NUMBER  OF  VECTORS  THAT  HAVE 
CONVERGED. 

M THE  NUMBER  OF  COLUMNS  OF  Q TO  ITERATE  WITH. 

EPS  A CONVERGENCE  CRITERION. 

MAX IT  AN  UPPER  BOUND  ON  THE  NUMBER  OF  ITERATIONS 

THE  PROGRAM  IS  TO  EXECUTE. 

START  AN  INITIALIZING  SIGNAL.  IF  START  .LT.  0. 

0  IS  INITIALIZED  BY  ORTHOGONALIZING  A SET  OF 
RANDOM  VECTORS.  IF  START  .GE.  0 THE  COLUMNS 
OF  Q ARE  USED  AS  A STARTING  APPROXIMATION  AND 
IF  START  .GE.  1 THEY  ARE  ALSO  ASSUMED  TO  BE 
ORTHONORMAL . 

*T  ON  RETURN  T CONTAINS  THE  REPRESENTATION  OF  A 

DESCRIBED  ABOVE. 

*ER  AN  ARRAY  THAT  ON  RETURN  CONTAINS  THE  REAL  PARTS 

OF  THE  EIGENVALUES  OF  T. 

BE I AN  ARRAY  THAT  ON  RETURN  CONTAINS  THE  COMPLEX  PARTS 

OF  THE  EIGENVALUES  OF  T. 

♦TYPE  AN  INTEGER  ARRAY.  ON  RETURN  TYPE(L)  CONTAINS 

0 IF  THE  L-TH  EIGENVALUE  IS  REAL 

1 IF  THE  L-TH  EIGENVALUE  IS  THE  FIRST  OF 
A COMPLEX  CONJUGATE  PAIR. 

2 IF  THE  L-TH  EIGENVALUE  IS  THE  SECOND  OF 
A COMPLEX  CONJUGATE  PAIR. 

-1  IF  THE  L-TH  EIGENVALUE  WAS  NOT  CORRECTLY 
DETERMINED. 

*RSD  AN  ARRAY  THAT  ON  RETURN  CONTAINS  THE  2-NORMS  OF 

THE  RESIDUAL  VECTORS  ABQ(BrL)  - Q*T<*»L>. 

BRSDX  AN  INTEGER  ARRAY  THAT  ON  RETURN  CONTAINS 

THE  ITERATIONS  AT  WHICH  THE  RESIDUALS  WERE  COMPUTED. 

WRITE  A LOGICAL  PARAMETER  THAT t IF  TRUE.  CAUSES 

INFORMATION  ABOUT  THE  COURSE  OF  THE  ITERATION  TO  BE 
WRITED  ON  UNIT  6. 


CONTROL  PARAMETERS 
INTEGER  INIT .SEED 

REAL  ALPHA . BET A . CNVTOL . GRPTOL . ORTTOL . STPF AC 
INTERNAL  VARIABLES 

I NTEGER  DORT . DSRR . I . I T . J . L . NGRP . NOGRP . NXTORT . 

1 NXTSRR .ORSDX ( 10 ) » OTYPE ( 10  > 

REAL  AE . AOE . AORSD » ARSD.CTR.OCTR.OEI ( 10) . OER < 1 0 > . ORSD ( 1 0 ) 

INITIALIZE  CONTROL  PARAMETERS 

INIT  - S 
STPF AC  - 2. 

SEED  - 69 
ALPHA  - 1. 

BETA  -1.1 
GRPTOL  - .001 
CNVTOL  - .001 
ORTTOL  - 2. 

INITIALIZE 


L 


1 
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IT  * 0 
DO  10  J=lrM 
RSD(J>  = 0. 

RSDX(J)  = -1 
TYPE  < J ) = -1 
10  CONTINUE 

IF< START  .GE.  0)  GO  TO  40 
DO  30  J=1»M 
DO  20  1=1 ,N 

Q ( I » J ) = RANDOM (SEED  > 

20  CONTINUE 
30  CONTINUE 
40  CONTINUE 

IF (START  .GT.  0)  GO  TO  50 
CALL  ORTH( Q» 1 > M* N ) 

50  CONTINUE 

SRR  LOOP 

100  CONTINUE 

IF  (WRITE)  WRITE ( 6 > 2000 ) IT,L 
2000  FORMAT (/10H  SRR  IT  =»I5»5H  L = .I3> 

CALL  SRRSTP (Of AO fATQ.L.MrN.T.ER.EI. TYPE t OER , OEI . OT YPE . 
1 WRITE) 

CALL  RESID( Q. AQ.  T . RSD.RSDX  »ORSD  »ORSDX  .L.M.M.N.IT.  TYPE  r 
1 WRITE) 


TEST  FOR  CONVERGENCE 
110  CONTINUE 

CALL  GROUP ( ER  > El r TYPE  t RSD  > GRPTOL  »L  »M  * N r 
1 NGRP  f CTR i AE  t ARSD» WRITE  > 

CALL  GROUP ( OER » OEI t OTYPE j ORSD » GRPTOL t L » M » N » 

1 NOGRP » OCTR  f AOE » AORSD  t WRITE ) 

IF ( NGRP  .NE.  NOGRP)  GO  TO  130 
IF(NGRP  .EQ.  0)  GO  TO  130 

IF ( ABS < AE-AOE ) .GT.  CTR*CNVTOL*FLOAT ( RSDX ( L ) -ORSDX ( L ) > ) 

1 GO  TO  130 

IF ( ARSD  .GT.  EPS)  GO  TO  130 
L = L + NGRP 
IF(L  .GT.  M)  GO  TO  130 
GO  TO  110 
130  CONTINUE 

IF  (WRITE)  WRITE ( 6 » 2000 ) IT.L 

EXIT  IF  THE  REQUIRED  NUMBER  OF  VECTORS  HAVE  CONVERGED. 

IF ( L .GT.  NV > GO  TO  300 

EXIT  IF  ITERATION  COUNT  EXCEEDS  THE  MAXIMUM  NUMBER 
OF  ITERATIONS. 

IF( IT  .GE.  MAXIT)  GO  TO  300 

DETERMINE  WHEN  THE  NEXT  SRR  STEP  IS  TO  BE  TAKEN. 

NXTSRR  = AMAX1  ( STPFAOFLOAT  ( IT ) » FLOAT  ( INIT ) ) 

NXTSRR  = M I NO ( MAXIT .NXTSRR ) 

DSRR  = NXTSRR-IT 
IF (NGRP  .NE.  NOGRP)  GO  TO  150 
IF (NGRP  .EO.  0)  GO  TO  150 
IF ( ARSD.GE . AORSD ) GO  TO  150 

DSRR  « ALPHA  + BETABFLOAT ( ORSDX < L ) -RSDX ( L )> BALOG ( ARSD/EPS > / 
1 ALOG( ARSD/ AORSD) 

DSRR  = MAXO(l.DSRR) 

150  CONTINUE 
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NXTSRR  » M I NO < NXTSRR  f IT+DSRR ) 

DETERMINE  THE  INTERVAL  BETWEEN  ORTHOGONAL I Z AT IONS 

DORT  - AMAX1 ( 1 . f0RTT0L/AL0G10(C0ND( T fMfWRITE) ) > 

NXTORT  » MINO(IT+DORTf NXTSRR) 

IF  (WRITE)  WRITE ( 6f  2001 ) IT fNXTSRR .NXTORT 
001  FORMAT < /10H  SRR  IT  »fI5f10H  NXTSRR  =fI5f10H  NXTORT  *fI3> 
DO  137  J*LfM 
DO  153  I-l.N 

Q( I r J)  = AO(IfJ) 

133  CONTINUE 
137  CONTINUE 
IT  - IT+1 

ORTHOGONAL I ZAT I ON  LOOP. 

160  CONTINUE 

POWER  LOOP 

170  CONTINUE 

IF< IT  ,E0.  NXTORT)  GO  TO  200 
CALL  ATQ<  Qf  AQ  fL  fM ) 

DO  190  J=LfM 
DO  180  1*1  fN 

Q(IfJ)  = AO(IfJ) 

180  CONTINUE 

190  CONTINUE 

IT  * IT  + 1 
GO  TO  170 
200  CONTINUE 

CALL  ORTH(QfLfMfN) 

NXTORT  = MINOdT+DORT  .NXTSRR) 

IF < IT  .LT.  NXTSRR)  GO  TO  160 
GO  TO  100 
300  CONTINUE 
NV  = L-l 
RETURN 
END 


SUBROUTINE  SRRSTP<  Of  AQf  ATQ  fL  fM  f N f T f ER  f El  f TYPE  f 
1 OERfOEIfOTYPEfWRITE) 

PARAMETERS  IN  THE  CALLING  SEQUENCE. 

INTEGER  LfMfNfTYPE<10)f0TYPE<10) 

REAL  Q(500f10) fAQ(500f 10) fT(10f 10) fER<10) fEI(IO) f 
1 OER(IO)fOEKIO) 

LOGICAL  WRITE 
EXTERNAL  ATQ 

SRRSTP  PERFORMS  A SCHUR-RA YLEIGH-RITZ  REFINEMENT  ON 
THE  SET  OF  M ORTHONORMAL  N-VECTORS  CONTAINED  IN 
THE  ARRAY  G.  FIRST  THE  SUBROUTINE  ATQ  IS  CALLED 
TO  GENERATE  THE  PRODUCT  OF  THE  MATRIX  A AND  THE 
VECTORS  Q IN  THE  ARRAY  Q.  THEN  THE  MATRIX 
T-TR(Q)*AQ  IS  REDUCED  TO  ORDERED  QUASI -TRIANGULAR 
FORM  BY  THE  SUBROUTINE  OTHESf  ORTRAN  AND  HQR3 . 

THE  REDUCING  TRANSFORMATION  V IS  POSTMULTIPLIED 
INTO  Q AND  AQ  TO  GIVE  THE  REFINDED  VECTORS  IN  Q AND 
THEIR  PRODUCT  WITH  A IN  AQ.  IT  IS  ASSUMENLD  THAT  IT  IS 
NECESSARY  TO  WORK  WITH  ONLY  COLUMNS  L THROUGH  M OF  T . 
THE  INFORMATION  CONTAINED  IN  POSITIONS  L THROUGH  M 
OF  THE  ARRAYS  ERf  EIf  AND  TYPE  IS  STORED  IN  THE 
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CORRESPONDING  POSITIONS  OF  THE  ARRAYS  CERf  OEI  AND  OTYPE. 
INTERNAL  VARIABLES 
INTEGER  I v J»K 

REAL  AP< IO)fP(IO) fVOOfIO) f MCHEPS 

IF  (WRITE)  WRITE ( 6 • 2000 ) L 
FORMAT (/12H  SRRSTP  L *fI3> 

MCHEPS  * 1. 

CONTINUE 

IF < MCHEPS+1 • .EQ.  I.)  GO  TO  5 
MCHEPS  » MCHEPS/2. 

GO  TO  3 
CONTINUE 

SAVE  THE  OLD  EIGENVALUES. 

DO  10  J»LfM 

OER ( J > = ER(J> 

OEKJ)  » EI(J> 

OTYPE(J)  - TYPE ( J > 

CONTINUE 

CALCULATE  THE  NEW  T. 

CALL  ATQ( 0* AQfL fM) 

DO  40  J»LfM 
DO  30  1*1  fM 
T(IfJ)  * 0. 

DO  20  K-1fN 

T(IfJ)  » T(IfJ)  + 0(KfI)*A0(KfJ) 

CONTINUE 

CONTINUE 

CONTINUE 

TRIANGULARIZE  T 


CALL  ORTHES(IOfMfLfMfTfP) 

CALL  ORTRAN(IOfMfLfMfTfPfV) 

CALL  HQR3< T f VfM fL f Mf MCHEPS f ERfEIf TYPE f 10 f 10) 
IF  (.NOT. WRITE)  GO  TO  48 
WRITE(6f1001) 

DO  43  1*1 fM 

WRITE ( 6 f 1000)  <T(IfJ)fJ=»1fM) 

CONTINUE 
WRITE (6f 1002 ) 

DO  43  1*1 fM 

WRITE ( 6 f 1000 ) (V(IfJ)fJ*1fM> 

CONTINUE 
WRITE(6f 1003) 

WRITE(AfIOOO)  (ER(I)fI-IfM) 

WRITE ( 6 f 1 004  > 

WRITE(AfIOOO)  (EI(I)fI-IfM) 

WRITE ( 6f 1005  > 

WRITE ( 6 f 1 100 ) ( TYPE ( I ) f 1*1 fM) 

FORMAT ( /1H  f10E12.4) 

FORMAT (//2H  T) 

FORMAT (//2H  V) 

FORMAT (//3H  ER) 

FORMAT (3H  El) 

FORMAT (//3H  TYPE) 

FORMAT ( /1H  f 101 12 ) 

CONTINUE 


r 
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TRANSFORM  Q AND  AQ 

DO  80  I-l.N 
DO  60  J»L»M 
P(J)  = 0. 

AP(J)  * 0. 

DO  50  K=1 f N 
P(J)  = P(J)  + 


Q(I.K)*V(K»J) 


AP(J>  + AQ(I,K)*V(K,J) 


C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


AP<  J) 

50  CONTINUE 
60  CONTINUE 

DO  70  J=L>M 
Q(I»J>  - P(J) 
AQ(I.J)  = AP(J) 
70  CONTINUE 
80  CONTINUE 
RETURN 
END 


SUBROUT  I NE  GROUP  (ERrEI > TYPE>RSD*  GRPTOL  iLtHiNf  NGRP  * CTR , AE » ARSD » 
I WRITE) 

PARAMETERS  IN  THE  CALLING  SEQUENCE. 

INTEGER  TYPE(IO)  »L»M»N»NGRP 

REAL  ERUO)  .El  (10)  ,RSDUO)>  GRPTOL.  CTR.  AEr  ARSD 
LOGICAL  WRITE 

GROUP  IS  A SUBROUTINE  TO  FIND  A CLUSTER  OF  COMPLEX 
NUMBERS  WHOSE  REAL  PARTS  ARE  CONTAINED  IN  THE  ARRAY 
ER  AND  IMAGINARY  PARTS  ARE  CONTAINED  IN  THE  ARRAY  El. 

THESE  NUMBERS  ARE  ASSUMED  TO  BE  STORED  IN  DESCENDING 
ORDER  OF  MAGNITUDE.  NGRP  IS  DETERMINED  AS  THE  LARGEST 
INTEGER  LESS  THAN  OR  EOUAL.  TO  M FOR  WHICH  THE  ABSOLUTE 
VALUE  E(J>  OF  THE  NUMBER  ER< J ) + El < J >*I  SATISFIES 


i 


E<L)  - E < L+NGRP-1 ) 


GRPTOL  / 


AND  FOR  WHICH  TYPE (L ) > TYPE (L+l )»...» TYPE ( L+NGRP-1 ) IS 
NONNEGATIVE.  IF  NGRP-O.  THE  SUBROUTINE  RETURNS 
CTR-AE-ARSD-O.  IF  NGRP .NE .0 » CTR  IS  SET  TO 
( E ( L ) +E ( L +NGRF - 1 ) ) /2 r AE  TO  THE  AVERAGE  OF  THE 
NUMBERS  ER+EI*I t AND  ARSD  TO  THE  RMS  AVERAGE  OF 
RSD<L>  fRSD<L+1 ) » . . . *RSD< L+NGRP-1 ) . 


INTERNAL  VARIABLES. 

INTEGER  JfLl 
REAL  MOD t MODI 
NGRP  - 0 

MOD  - SORT (ER(L  >**2  + EI(L)**2) 

CTR  » 0. 

10  CONTINUE 

LI  ■ L + NGRP 

IF(L1.GT.M  .OR.  TYPE<Ll).LT.O)  GO  TO  20 
MODI  » SORT <ER(L1 ) **2  + EI(L1)**2) 

IF <ABS< MOD-MODI)  .GT.  GRPTOL* < M0D+M0D1 ) ) GO  TO  20 
CTR  - (MOD  ♦ MODI >/2. 

NGRP  - NGRP  ♦ TYPE (LI ) + 1 
GO  TO  10 
20  CONTINUE 
AE  - 0. 

ARSD  - 0. 


L 
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IF < NGRP  .EQ.  0)  GO  TO  40 
LI  « L+NGRP-1 
DO  30  J»LfL1 

AE  =■  AE  + ER  < J ) 

ARSD  = ARSD  + RSD(J)**2 
30  CONTINUE 

AE  » AE/FLOAT  ( NGRP ) 

ARSD  » SORT (ARSD/FLOAT< NGRP ) ) 

40  CONTINUE 

WRITE ( 6 f 2000 > NGRP fCTRfAE f ARSD 
2000  FORMAT </14H  GROUP  NGRP  =fI3f7H  CTR  =fE12.4f6H  AE  »fE12.4f 
1 8H  ARSD  ” f £12 . 4 ) 

RETURN 

END 


I 


SUBROUTINE  RESID < Q f AO r T fRSDf RSDX fORSD f ORSDX f LI f L2 f M > N f IT f TYPE f 
1 WRITE) 

C 

C PARAMETERS  IN  THE  CALLING  SEQUENCE. 

C 

INTEGER  RSDX  (10)  fORSDX  ( 10  ) fLI  fL2  fMfNf  IT  f TYPE  < 10 ) 

REAL  Q<300f10>fAQ(S00f10)fT<10f10>fRSD<10)f0RSD(10> 

LOGICAL  WRITE 
C 

C RESID  COMPUTES  RESIDUALS  CORRESPONDING  TO  EIGENVALUES 

C LI  THROUGH  L2  OF  THE  QUASI -TRIANGULAR  MATRIX  T OF 
C ORDER  M.  SPECIFICALLY  IF  T(*fJ)  IS  THE  J-TH  COLUMN 

C OF  Tf  RSD(J)  IS  SET  TO  NORM ( AO <*fJ)-QT(*fJ))f  WHERE  THE 

C NORM  IS  THE  EUCLIDEAN  NORM.  THE  INDEX  RSDX(J)  IS  SET 

C EQUAL  TO  IT.  IF  THE  J-TH  EIGENVALUE  IS  COMPLEX 

C ( TYPE ( J ) - 1)  THE  RMS  AVERAGE  OF  RSD(J)  AND  RSD(J+1) 

C IS  PLACED  IN  RSD(J)  AND  RSD<J+1>.  THE  INITIAL  VALUES  OF 

C RSD  AND  RSDX  ARE  STORED  IN  ORSD  AND  ORSDX. 

C 

C 

C INTERNAL  VARIABLES 

C 

INTEGER  IfJfKfKU 
REAL  S 
C 

IF  (WRITE)  WRITE (6 f 2000 ) L1fL2 
2000  FORMAT ( /12H  RESID  LI  »fI5f6H  L2  «fI5> 

IF (LI  .GT.  L2)  RETURN 
DO  30  J-LIfI.2 

ORSD(J)  = RSD ( J > 

ORSDX ( J ) » RSDX(J) 

RSDX(J)  » IT 
KU  - MIN0(J+1fM> 

IF ( TYPE ( J ) .EQ.  0)  KU  = J 
RSD(J)  • 0. 

DO  20  I-1fN 
S » 0. 

DO  10  K“1 »KU 

S » S + Q(IfK)»T(KfJ) 

10  CONTINUE 

RSD(J)  - RSD(J)  * (AQ(IfJ)-S) **2 
20  CONTINUE 
30  CONTINUE 

DO  30  J-L1fL2 

IF( TYPE ( J)  .NE.  1)  GO  TO  40 

RSD(J>  « (RSD(J)  4 RSD ( J+l ) ) /2 . 

RSD ( J+ 1 ) - RSD(J) 

40  CONTINUE 

RSD(J)  - SORT (RSD( J) ) 

SO  CONTINUE 


i 
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IF  (.NOT.  WRITE)  GO  TO  60 
WRITE <6f 1000) 

WRITE(6f 1001 ) <RSD< I ) » 1*1 .M) 
WRITE (6f 1002) 

WRITE (6f 1003 ) < RSDX  <I)»I»1.M) 

1000  FORMAT <//4H  RSD ) 

1001  FORMAT </lH  f10E12.4> 

1002  FORMAT < //5H  RSDX) 

1003  FORMAT < /1H  rl0I12) 

C 

60  CONTINUE 
RETURN 
END 


C 

C 

C 


C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


100 

c 

c 

c 


110 


c 

c 

c 

c 

c 

c 


SUBROUTINE  ORTH(QfLfMfN) 

PARAMETERS  IN  THE  CALLING  SEQUENCE 

INTEGER  LfM.N 
REAL  Q<S00f 10) 

ORTH  ORTHONORMAL I ZES  COLUMNS  L THROUGH  M OF  THE  ARRAT 
Q WITH  RESPECT  TO  COLUMN  1 THROUGH  M.  COLUMNS  1 
THROUGH  L-l  ARE  ASSUMED  TO  BE  ORTHONORMAL.  THE  METHOD 
IS  THE  GRAM-SCHMIDT  METHOD  WITH  REORTHOGONCELIZATION. 

A COLUMN  IS  ACCEPTED  WHEN  AN  ORTHOGONAL I ZAT I ON  DOES 
NOT  REDUCE  ITS  EUCLIDEAN  NORM  BY  A FACTOR  OF  MORE 
THAN  TOL.  IF  THIS  IS  NOT  DONE  IN  MAXTRY  ATTEMPTS 
THE  PROGRAM  STOPS.  THE  PROGRAM  ALSO  STOPS  IF  IT 
ENCOUNTERS  A ZERO  VECTOR. 

INTERNAL  CONTROL  PARAMETERS 

REAL  TOL 
INTEGER  MAXTRY 

INTERNAL  VARIABLES 

REAL  NORM >00 

INTEGER  If JfJMIfKfTRY 

MAXTRY  - 5 

TOL  » .5 

DO  160  J*L»M 

ORTHOGONALIZE  THE  J-TH  VECTOR 


JM1  * J-l 
TRY  - 0 
CONTINUE 

COMPUTE  THE  NORM  OF  THE  VECTOR. 

NORM  - 0. 

DO  110  1*1 f N 

NORM  - NORM  4 Q<IfJ)**2 
CONTINUE 

NORM  ■ SORT (NORM) 

ERROR  TEST 

IF ( NORM  ,E0.  0.)  GO  TO  170 
SCALE  THE  VECTOR. 

00  120  1*1 »N 
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Q(IfJ)  = Q(IfJ)/NORM 
120  CONTINUE 

TEST  TO  SEE  IF  THE  J-TH  VECTOR  IS  ORTHOGONAL. 

IF(J  .EQ.  1)  GO  TO  160 
IF (TRY  .EQ.  0)  NORM  * 0. 

IF (NORM  .GT.  TOL)  GO  TO  160 
TRY  ■ TRY  + 1 

IF ( TRY  .GT.  MAXTRY)  GO  TO  170 

PERFORM  ONE  MODIFIED  GRAM-SCHMIDT  STEP. 

DO  ISO  K-1.JM1 
QQ  - 0. 

DO  130  I»1fN 

QQ  = QQ  + Q(I'K>*Q(IfJ> 

130  CONTINUE 

DO  140  I*1fN 

Q(ItJ)  = Q(IfJ)  - QQ*Q(IrK) 

140  CONTINUE 

150  CONTINUE 
GO  TO  100 
160  CONTINUE 
RETURN 

170  CONTINUE 

WRITE ( 6 f 2000  > 

000  FORMAT ( /I 4H  ERROR  IN  ORTH) 

STOP 

END 


REAL  FUNCTION  COND( TfMfWRITE) 

PARAMETERS  IN  THE  CALLING  SEQUENCE 

REAL  T( 10f 10 ) 

INTEGER  M 
LOGICAL  WRITE 

CCND  IS  A FUNCTION  THAT  RETURNS  THE  CONDITION 
NUMBER  WITH  RESPECT  TO  THE  ROW-SUM  NORM  OF  THE  UPPER 
HESSENBERG  MATRIX  T OF  ORDER  M. 


INTERNAL  VARIABLES 

REAL  MULT ( 10) fNT fNTRfNTI fNTIRfTKIOfIO) 

INTEGER  IfIIfJfJMIfJIfKfPVT(IO) 

MM1  » M-l 
NT  » 0. 

DO  20  I«1fM 

II  » MAXO(I-IfI) 

NTR  » 0. 

DO  10  J-UfM 

TI(IfJ)  = T( I f J> 

NTR  » NTR  + ABS(T(IfJ>) 

10  CONTINUE  \ 

NT  - AMAXKNTfNTR) 

20  CONTINUE 

DO  60  I«1fMM1 
PVT ( I ) » 0 
MULT(I)  » 0. 

IF(T1(I+1fI)  .EQ.  0.)  GO  TO  60 

IF( ABS( T1 (I+1fI)). LE . ABS(T1(IfI>)>  GO  TO  40 


PVT  < I ) - 1 
DO  30  J-IfM 
S ■ T1 ( I • J) 

T1(IfJ>  « Tl(ltl,J) 

TKI  + l.J)  » S 
30  CONTINUE 
40  CONTINUE 

MULT(l)  - TKItl.D/TKI.I) 

T1 < 1 + 1 » l ) » 0. 

II  » 1 + 1 

DO  30  J-I1fH 

TKI  + l.J)  » THI  + l.J)  - MULT(I)*T1(IfJ) 

30  CONTINUE 
60  CONTINUE 

DO  110  J»1  fM 

IF(T1(JfJ)  .NE.  0.)  GO  TO  70 
COND  = 1.E8 
RETURN 
70  CONTINUE 

Tl(J.J)  =■  1./T1(JfJ> 

IF < J .EQ.  1)  GO  TO  100 
JH1  = J-l 
DO  90  I-1.JM1 
S * 0. 

DO  80  K=IfJM1 

S * S + T1(IfK)*T1(KfJ> 

80  CONTINUE 

T1(IfJ>  ■ -SITl(JiJ) 

90  CONTINUE 

100  CONTINUE 
110  CONTINUE 

DO  160  JJ»1fMM1 
J » M-JJ 
J1  » J+l 

IF ( MULT ( J ) .EQ.  0.)  GO  TO  130 
DO  120  I-1fJ1 

Tl(IiJ)  » TI(IfJ)  - MULT(J)*T1(IfJ+1) 

120  CONTINUE 

130  CONTINUE 

IF < PVT < J ) .EO.  0)  GO  TO  150 
DO  140  I-1.J1 
S =»  Tl(liJ) 

TI(IfJ)  » T1 < I f J+l > 

T1(IfJ+1)  » S 
140  CONTINUE 

130  CONTINUE 
160  CONTINUE 
NT1  = 0. 

DO  180  I»1fM 

INI  » MAXO(IfI-I) 

NT1R  - 0. 

DO  170  J-IM1fM 

NT1R  » NT 1R  ♦ ABS(T1(IfJ)> 

170  CONTINUE 

NT1  » AMAXKNTIfNTIR) 

180  CONTINUE 

COND  ■ NT4NT1 

IF  (WRITE)  WRITE ( 6 f 2000 ) NTfNTIfCOND 
2000  FORMAT < /11H  COND  NT  »fE12.4f6H  NT1» f E12 . 4 f 8H  COND  =*fE12.4> 
RETURN 
END 


FUNCTION  RANDOM (SEED) 

INTEGER  SEED 

C RANDOM  IS  A FUNCTION  THAT  PRODUCES  A PSEUDO-RANDOM 


r 
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FLOATING  POINT  NUMBER  IN  THE  INTERVAL  FROM  ZERO  TO  0!*E. 
SEED  - MOD  < 4621 tSEED+21 1 3 > lOOOO ) 

RANDOM  » FLOAT ( SEED ) /I .E4 

RETURN 

END 


SUBROUTINE  HQR3( A> V ,N , NLOH.NUP f EPS .ER , El , TYPE , NA.NV) 

C 

INTEGER  N » NA , NLOU » NUP  » NV » TYPE  < N > 

REAL  A < NA ,N),EI(N) >ER <N ) . EPS » V ( NV , N ) 

C 

C HQR3  REDUCES  THE  UPPER  HESSENBERG  MATRIX  A TO  QUASI- 
C TRIANGULAR  FORM  BY  UNITARY  SIMILARITY  TRANSFORMATIONS. 

C THE  EIGENVALUES  OF  A,  WHICH  ARE  CONTAINED  IN  THE  1X1 

C AND  2X2  DIAGONAL  BLOCKS  OF  THE  REDUCED  MATRIX.  ARE 

C ORDERED  IN  DESCENDING  ORDER  OF  MAGNITUDE  ALONG  THE 

C DIAGONAL.  THE  TRANSFORMATIONS  ARE  ACCUMULATED  IN  THE 

C ARRAY  V.  H0R3  REQUIRES  THE  SUBROUTINES  EXCHNG. 

C QRSTEP.  AND  SPLIT.  THE  PARAMETERS  IN  THE  CALLING 

C SEQUENCE  ARE  (STARRED  PARAMETERS  ARE  ALTERED  BY  THE 

C SUBROUTINE) 

C 

C *A 

C 
C 

c 

c tv 

c 
c 

C N 

C NLOW 

C NUP 

c 
c 

C 

c 

C EPS 

C tER 

c 

C tEI 

C 

C tTYPE 

C 

c 
c 
c 
c 
c 
c 

C NA 

C NV 

c 

C INTERNAL  VARIABLES 

C 

INTEGER  I.IT.L.MU.NL.NU 
REAL  E1.E2.P.Q.R.S.T.W.X.Y.Z 
LOGICAL  FAIL 
C 

C INITIALIZE. 

C 

DO  10  I -NLOW. NUP 
TYPE<I)  - -1 
10  CONTINUE 
T - 0. 

C 


AN  ARRAY  THAT  INITIALLY  CONTAINS  THE  N X N 
UPPER  HESSENBERG  MATRIX  TO  BE  REDUCED.  ON 
RETURN  A CONTAINS  THE  REDUCED.  QUASI- 
TRI ANGULAR  MATRIX. 

AN  ARRAY  THAT  CONTAINS  A MATRIX  INTO  WHICH 
THE  REDUCING  TRANSFORMATIONS  ARE  TO  BE 
MULTIPLIED. 

THE  ORDER  OF  THE  MATRICES  A AND  V. 

A ( NLOW- 1 .NLOW)  AND  A(NUF'.NUP+U>  ARE 
ASSUMED  TO  BE  ZERO.  AND  ONLY  ROWS  NLOW 
THROUGH  NUP  AND  COLUMNS  NLOU  THROUGH 
NUP  ARE  TRANSFORMED.  RESULTING  IN  THE 
CALCULATION  OF  EIGENVALUES  NLOU 
THROUGH  NUP. 

A CONVERGENCE  CRITERION. 

AN  ARRAY  THAT  ON  RETURN  CONTAINS  THE  REAL 
PARTS  OF  THE  EIGENVALUES. 

AN  ARRAY  THAT  ON  RETURN  CONTAINS  THE 
IMAGINARY  PARTS  OF  THE  EIGENVALUES. 

AND  INTEGER  ARRAY  WHOSE  I-TH  ENTRY  IS 

0 IF  THE  I-TH  EIGENVALUE  IS  REAL, 

1 IF  THE  I-TH  EIGENVALUE  IS  COMPLEX 
WITH  POSITIVE  IMAGINARY  PART. 

2 IF  THE  I-TH  EIGENVALUE  IS  COMPLEX 
WITH  NEGATIVE  IMAGINARY  PART, 

-1  IF  THE  I-TH  EIGENVALUE  WAS  NOT 
CALCULATED  SUCCESSFULLY. 

THE  FIRST  DIMENSION  OF  THE  ARRAY  A. 

THE  FIRST  DIMENSION  OF  THE  ARRAY  V. 
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C MAIN  LOOP.  PINO  AND  ORDER  EIGENVALUES. 

C 

NU  - NUP 

100  IF < NU  .LT.  NLOW > GO  TO  500 
IT  - 0 
C 

C OR  LOOP.  FIND  NEGLIGABLE  ELEMENTS  AND  PERFORM 

C QR  STEPS. 

C 

110  CONTINUE 

c 

C SEARCH  BACK  FOR  NEGLIGABLE  ELEMENTS. 

C 

L - NU 


C 

C 

C 

c 


c 

c 

c 

c 


c 

c 

c 


I 


c 

c 

c 

c 


1 

I 

\ 

\ 


120  CONTINUE 

IF<L  .EQ.  NLOW)  GO  TO  130 

IF(ABS<A<LfL-1> > .LT.  EPS«< ABS< A<L-1 fL-1 ) >4AB8<A<L'L) ) > ) 
1 GO  TO  130 

L » L-l 
GO  TO  120 
130  CONTINUE 

TEST  TO  SEE  IF  AN  EIGENVALUE  OR  A 2X2  BLOCK 
HAS  BEEN  FOUND. 

X = A(NU.NU) 

IF <L  .EQ.  NU)  GO  TO  300 
Y » A(NU-l.NU-l) 

W - A<NU.NU-1)*A<NU-1.NU) 

IF(L  .EQ.  NU-1)  GO  TO  200 

TEST  ITERATION  COUNT.  IF  IT  IS  30  QUIT.  IF 
IT  IS  10  OR  20  SET  UP  AN  AD-HOC  SHIFT, 

IF(  IT  .EQ.'  30)  GO  TO  500 
IF(IT,N£.10  .AND.  IT.NE.20)  GO  TO  150 

AD-HOC  SHIFT. 

T =*  T 4 X 
DO  14b  1»NL0UfNU 

A(IfI)  ■'  A < I . I ) - X 
140  ' CONTINUE 

S - ABS ( A < NU t NU- 1 > > 4 ABS< A<NU-1 fNU-2) ) 

X » 0.75*S 
Y « X 

W - -0.4375*S**2 
ISO  CONTINUE 

IT  - IT  + 1 

LOOK  FOR  TWO  CONSECUTIVE  SMALL  SUB-DIAGONAL 
ELEMENTS. 

NL»  NU-2 
160  CONTINUE 

Z - A<NLfNL) 

R - X - Z 

S * Y - Z 

P » ( R4S-W) /A CNL  + 1 »NL ) 4 A<NLfNL41) 

Q - A<NL4l\NL41>  - z - r - s 
R ■ A<  NL+2  fNL  + 1 ) 

S - ABS(P)  4 ABS (Q > 4 ABS(R) 

P - P/S 

a - a/s 

R - R/S 

IF<NL  .EQ. 


L>  GO  TO  170 


r 


y. 


r s 


i 


l! 


! 


! 


V 
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IF(ABS(A<NLfNL-1> ) * < ABS < Q > +ABS < R > ) .LE. 

1 EPS*ABS<P)*<ABS< A<NL-1 fNL-1 > > +ABS < Z >+ABS < A< NL  + 1 f NL  + 1 ) ) > ) 

2 GO  TO  170 
NL  - NL-1 

GO  TO  160 
170  CONTINUE 
C 

C PERFORM  A OR  STEP  BETWEEN  NL  AND  NU. 

C 

CALL  QRSTEP<AfVfPfQfRfNLfNUfNfNAfNV> 

GO  TO  110 

c 

C 2X2  BLOCK  FOUND. 

C 

200  IF (NU  .NE.  NLOW+1 ) A<NU-1fNU-2>  = 0. 

A(NUfNU)  = A(NU.NU)  + T 
A(NU-IfNU-I)  = A ( NU-1 f NU— 1 ) + T 
TYPE ( NU)  = 0 
TYPE < NU-1 ) * 0 
MU  = NU 
C 

C LOOP  TO  POSITION  2X2  BLOCK. 

C 

210  CONTINUE 

NL  » MU-1 
C 

C ATTEMPT  TO  SPLIT  THE  BLOCK  INTO  TWO  REAL 

C EIGENVALUES. 

C 

CALL  SPLIT <AfVfNfNLfE1fE2fNAfNV) 

C 

C IF  THE  SPLIT  WAS  SUCCESSFUL r GO  AND  ORDER  THE 

C REAL  EIGENVALUES. 

C 

IF<  A<  MUfMU-1 ) .Ed.  0.)  GO  TO  310 
C 

C TEST  TO  SEE  IF  THE  BLOCK  IS  PROPERLY  POSITIONED. 

C AND  IF  NOT  EXCHANGE  IT 

C 

IF < MU  .EQ.  NUP)  GO  TO  400 
IF < MU  .EQ.  NUP-l)  GO  TO  220 
IF<  A<  HU+2  fMU+1 ) .EQ.  0.)  GO  TO  220 
C 

C THE  NEXT  BLOCK  IS  2X2. 

C 

IF  < A < MU—  1 1 MU- 1 ) #A  ( MU f MU ) -A  ( MU- 1 f MU ) *A  ( MU f MU—  1 ) 

1 .GE.  A<MU+1fMU+1)*A<MU+2fMU+2)-A<MU+1fMU+2>* 

2 A < HU+2 fMU+1 ) ) 

3 GO  TO  400 

CALL  EXCHNG<AfVfNfNLf2f2fEPSfFA1LfNAfNV> 

IF < .NOT . FAIL)  GO  TO  215 
TYPE  < NL ) - -1 
TYPE < NL+1 ) » -1 
TYPE(NL+2)  * -1 
TYPE<NL+3)  » -1 
GO  TO  500 

215  CONTINUE 

MU  - MU+2 

GO  TO  230 

220  CONTINUE 
C 

C THE  NEXT  BLOCK  IS  1X1. 

C 

IF<A<MU-1fHU-1)»A<MUfMU>-A<MU-1fMU)*A<MUfNU-1> 

1 . GE.  A<MU+1 »MU+1 )*»2> 

2 GO  TO  400 
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c 

c 

c 


c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 


CALL  EXCHNG( Ar Vr Nr NL  r 2r  1 r EPS r FAIL  r NArNV ) 

IF ( . NOT • FAIL)  GO  TO  223 
TYPE(NL)  » -1 
TYPE ( NL+1 ) « -I 
TYPE<NL+2 ) » -1 
GO  TO  500 

223  CONTINUE 

NU  « HU+1 
230  CONTINUE 
GO  TO  210 

SINGLE  EIGENVALUE  FOUND. 

300  NL  » 0 

A < NU  r NU ) - A < NU  r NU  > + T 

IF( NU  .NE.  NLOU)  A(NUrNU-l)  = 0. 

TYPE(NU)  » 0 
MU  » NU 

LOOP  TO  POSITION  ONE  OR  TWO  REAL  EIGENVALUES. 

310  CONTINUE 

POSITION  THE  EIGENVALUE  LOCATED  AT  A(NLrNL) • 

320  CONTINUE 

IF  < MU  .EQ.  NUP)  GO  TO  350 

IF (MU  .EQ.  NUP-1)  GO  TO  330 

IF ( A(MU+2 r MU+1  ) .EQ.  0.)  GO  TO  330 

THE  NEXT  BLOCK  IS  2X2. 

IF(A(MUrMU>**2  .GE. 

1 A(MU+lrMU+l )*A(MU+2rNU+2)-A<NU+lrMU+2>*A(MU+2rMU+l) ) 

2 GO  TO  400 

CALL  EXCHNG( At VrNrMUr  1 r 2rEPSrFAILrNArNV) 

IF(.NOT.  FAIL)  GO  TO  323 
TYPE (MU)  ■ -1 
TYPE (MU+1 ) » -1 
TYPE ( MU+2 ) = -1 
GO  TO  500 

323  CONTINUE 

MU  - MU+2 
GO  TO  340 
330  CONTINUE 

THE  NEXT  BLOCK  IS  1X1. 

IF(ABS(A(MUrMU) > .GE.  ABS ( A(MU+1 r MU+1 ) ) > 

X GO  TO  350 

CALL  £XCHNG( Ar VrNrMUr 1 » 1 rEPSrFAILrNArNV) 

MU  » MU+1 

340  CONTINUE  ' 

GO  TO  320 
350  CONTINUE 

MU  - NL 
NL  - 0 

IF ( MU  .NE.  0)  GO  TO  310 

GO  BACK  AND  GET  THE  NEXT  EIGENVALUE. 

400  CONTINUE 
NU  - L-l 
GO  TO  100 

ALL  THE  EIONVALUES  HAVE  BEEN  FOUND  AND  ORDERED. 


r 
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COMPUTE  THEIR  VALUES  AND  TYPE. 

500  IF ( NU  .LT.  NLOU)  GO  TO  507 
DO  503  1 3 1 1 NU 

A(I»I>  = Adfl)  + T 
503  CONTINUE 
507  CONTINUE 
NU  » NUP 
510  CONTINUE 

IF ( TYPE ( NU ) .NE.  -1)  GO  TO  515 
NU  = NU-1 
GO  TO  540 
515  CONTINUE 

IF ( NU  .EQ.  NLOU)  GO  TO  520 
IF(A(NU»NU-1 > . EQ.  0.)  GO  TO  520 

2X2  BLOCK. 

CALL  SPLIT ( A» V»Nf NU-1 .El .E2»NA»NV> 
IF  ( A ( NU f NU- 1 ) .EQ.  0.)  GO  TO  520 
ER<NU)  = El 
El (NU-1)  = E2 
ER(NU-l)  = ER(NU) 

E I < NU ) = -EI(NU-l) 

TYPE(NU-l)  - 1 
TYPE(NU)  = 2 
NU  = NU-2 
GO  TO  530 
520  CONTINUE 

SINGLE  ROOT. 

ER(NU)  = A(NUrNU) 

E I < NU ) » 0. 

NU  - NU-1 
530  CONTINUE 
540  CONTINUE 

IF ( NU  .GE.  NLOU)  GO  TO  510 

RETURN 

END 
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SUBROUTINE  EXCHNG ( Ar V, N. L , B1 . B2» EPS. FAIL r NA ,NV> 

INTEGER  81 » B2 » L . NA . NV 
REAL  A(NA.N) .EPS»V(NV.N) 

LOGICAL  FAIL 

GIVEN  THE  UPPER  HESSENBERG  MATRIX  A UITH  CONSECUTIVE 
B1XB1  AND  B2XB2  DIAGONAL  BLOCKS  (B  . B2  ,3E.  2) 
STARTING  AT  A(LrL) r EXCHNG  PRODUCES  A UNITARY 
SIMILARITY  TRANSFORMATION  THAT  EXCH+5GES  THE  B36B2S 
ALONG  WITH  THEIR  EIGENVALUES.  THE  TRANSFORMATION 
IS  ACCUMULATED  IN  V.  EXCHNG  REQUIRES  THE  SUBROUTINE 
QRSTEP.  THE  PARAMETERS  IN  THE  CALLING  SEQUENCE  ARE 
(STARRED  PARAMETERS  ARE  ALTERED  BY  THE  SUBROUTINE) 

*A  THE  MATRIX  UHOSE  BLOCKS  ARE  TO  BE 

INTERCHANGED. 

*V  THE  ARRAY  INTO  UHICH  THE  TRANSFORMATIONS 

ARE  TO  BE  ACCUMULATED. 

N THE  ORDER  OF  THE  MATRIX  A. 

L THE  POSITION  OF  THE  BLOCKS. 

B1  THE  SIZE  OF  THE  FIRST  BLOCK. 

B2  THE  SIZE  OF  THE  SECOND  BLOCK. 

EPS  A CONVERGENCE  CRITERION. 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 

c 

c 


3 


3 


7 


c 

c 

c 


10 


20 


30 


•FAIL  A LOGICAL  VARIABLE  WHICH  IS  FALSE  ON  A 

NORMAL  RETURN.  IF  THIRTY  ITERATIONS  WERE 
PERFORMED  WITHOUT  CONVERGENCE f FAIL  IS  SET 
TO  TRUE  AND  THE  ELEMENT 
A<L+B2fL+B2-1 > CANNOT  BE  ASSUMED  ZERO. 

HA  THE  FIRST  DIMENSION  OF  THE  ARRAY  A. 

NV  THE  FIRST  DIMENSION  OF  THE  ARRAY  V. 

INTERNAL  VARIABLES. 

INTEGER  I f IT » J«L1 »M 
REAL  PfQfRfSfWfXfYfZ 

FAIL  = .FALSE. 

IFCBl  .EQ.  2)  GO  TO  40 
IF (B2  .EQ.  2)  GO  TO  10 

INTERCHANGE  1X1  AND  1X1  BLOCKS. 

LI  * L+l 

Q = A<L+1 »L+1 ) - A ( L fL  > 

P » A < L fL+1  ) 

R = AMAXKPfQ) 

IF (R  .EQ.  0.)  RETURN 
P = P/R 
Q » Q/R 

R = SORT ( P**2  + 0**2) 

P » P/R 
Q * Q/R 
DO  3 J-L fN 

S ■ P*A(LpJ)  + Q*A<L+1fJ> 

A(L+1fJ)  = P*A(L+1 » J)  - Q*A<Lf J) 

A ( L f J ) = S 
CONTINUE 
DO  3 1*1 fLI 

S - P*A< I fL)  + Q*A<IfL+1> 

A< I fL+1 ) = P*A<IfL+1)  - Q*A(IfL) 

A(IfL)  * s 
CONTINUE 
DO  7 1*1 fN 

S - P*V< I fL)  + 0*V(IfL+1) 

V<IfL+1)  = P*V(IfL+1)  - Q*V< IfL) 

V< I fL)  = S 
CONTINUE 
A(L+1fL)  * 0. 

RETURN 

CONTINUE 

INTERCHANGE  1X1  AND  2X2  BLOCKS. 

X - A(LfL) 

P - 1. 

Q * 1. 

R - 1. 

CALL  QRSTEP <AfVfPfQfRfLfL+2fNfNAf NV ) 

IT  - 0 
IT  - IT+1 

IF ( IT  ,LE.  30)  GO  TO  30 
FAIL  - .TRUE. 

RETURN 
CONTINUE 
P » A<LfL)  - X 
Q - A(L+1fL> 

R * 0. 

CALL  QRSTEP <AfVfPf0fRfLfL+2fNf NA f NV ) 

IF  < ABS(A(L+2fL+1 ) ) .GT. 
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1 EPS* ( ABS  (A(L+1»L+1) )+ABS  ( A (L+2»L+2) ) ) ) 

1 GO  TO  20 

A(L+2,L+1)  = 0, 

RETURN 
40  CONTINUE 
C 

C INTERCHANGE  2X2  AND  B2XB2  BLOCKS. 

C 

M = L+2 

IF ( B2  .EQ.  2)  M = M+l 
X = A<  L + l rL+1 ) 

Y - A<L»L> 

W = A(l.+  1.L)*A(L,L+1) 

P - 1. 

Q ■ 1. 

R » 1 . 

CALL  QRSTEP(A.VrP»Q,R,LrM»N»NA,NV> 

IT  = 0 

50  IT  = IT+1 

IF( IT  . LE.  30)  GO  TO  60 
FAIL  = .TRUE. 

RETURN 
60  CONTINUE 

Z = A(L»L) 

R = X - Z 
S » Y - Z 

P » (R*S-W)/A(L+1 *L)  + A<L»L+1 ) 

Q » A<L+1 »L+1 ) - Z - R - S 
R = A(L+2*L+1> 

S * ABS<P)  + ABS(Q)  + ABS(R> 

P - P/S 

0 = 0/S  \ 

R - R/S 

CALL  QRSTEPCA'V'P'Q.RrL’M'N'NArNV) 

IF< ABS< A<M-1 f N-2) ) .GT.  EPS* < ABS< A<M-1 , M-l ) ) +ABS< A(M-2»M-2) ) > ) 

1 GO  TO  50 

A<N-lfN-2)  = 0. 

RETURN 

CONTINUE 

END 


* 

■' 
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SUBROUTINE  SPLIT(ArV»NrL»El r E2.NA.NV) 

INTEGER  L.N'NA.NV 
REAL  A < NA  » N ) * U < NV  r N ) 

GIVEN  THE  UPPER  HESSENBERG  MATRIX  A WITH  A 2X2  BLOCK 
STARTING  AT  A(L»L) » SPLIT  DETERMINES  IF  THE 
CORRESPONDING  EIGENVALUES  ARE  REAL  OR  COMPLEX.  IF  THEY 
ARE  REALr  A ROTATION  IS  DETERMINED  THAT  REDUCES  THE 
BLOCK  TO  UPPER  TRIANGULAR  FORM  WITH  THE  EIGENVALUE 
OF  LARGEST  ABSOLUTE  VALUE  APPEARING  FIRST.  THE 
ROTATION  IS  ACCUMULATED  IN  V.  THE  EIGENVALUES  (REAL 
OR  COMPLEX)  ARE  RETURNED  IN  El  AND  E2.  THE  PARAMETERS 
IN  THE  CALLING  SEQUENCE  ARE  (STARRED  PARAMETERS  ARE 
ALTERED  BY  THE  SUBROUTINE) 

*A  THE  UPPER  HESSENVERG  MATRIX  WHOSE  2X2 

BLOCK  IS  TO  BE  SPLIT. 

*V  THE  ARRAY  IN  WHICH  THE  SPLITTING  TRANS- 

FORMATION IS  TO  BE  ACCUMULATED. 

N THE  ORDER  OF  THE  MATRIX  A. 

L THE  POSITION  OF  THE  2X2  BLOCK. 

*E1  ON  RETURN  IF  THE  EIGENVALUES  ARE  COMPLEX 

*E2  El  CONTAINS  THEIR  COMMON  REAL  PART  AND 
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E2  CONTAINS  THE  POSITIVE  IMAGINARY  PART. 
IF  THE  EIGENVALUES  ARE  REAL*  El  CONTAINS 
THE  ONE  LARGEST  IN  ABSOLUTE  VALUE  AND  E2 
CONTAINS  THE  OTHER  ONE. 

NA  THE  FIRST  DIMENSION  OF  THE  ARRAY  A. 

NV  THE  FIRST  DIMENSION  OF  THE  ARRAY  V. 

INTERNAL  VARIABLES 

INTEGER  I,J»L1 

REAL  PrQ.Rf TfU.UfXfYfZ 

X = A( L + l fL  + 1 ) 

Y - A<LfL) 

U = A<LfL+1 )*A<L+1 fL> 

P * <Y-X)/2. 

Q = P**2  + W 

IF(Q  .GE.  0. > GO  TO  5 

COMPLEX  EIGENVALUE. 

El  = P + X 
E2  = SORT  < -0 ) 

RETURN 
5 CONTINUE 

TUO  REAL  EIGENVALUES.  SET  UP  TRANSFORMATION. 

Z * SORT ( Q ) 

IF< P .LT.  0. ) GO  TO  10 
Z * P + Z 
GO  TO  20 
10  CONTINUE 

Z =•  P - Z 
20  CONTINUE 

IF < Z .EG.  0.)  GO  TO  30 
ft  * -W/Z 
GO  TO  40 
30  CONTINUE 
R » 0. 

40  CONTINUE 

IF<ABS(X+Z)  .GE.  ABS < X+R ) ) Z » R 

Y = Y - X - Z 
X * -Z 

T » A < L fL+1 ) 

U - A<L+lrL> 

IF(ABS< Y)+ABS(U)  .LE.  ABS  < T > +ABS < X ) > GO  TO  60 

0 » u 

P = Y 
GO  TO  70 
60  CONTINUE 
0 » X 
P =*  T 

70  CONTINUE 

R * SORT ( P**2  + 0**2) 

IF<R  .GT.  0.  ) GO  TO  80 
El  - A(LrL) 

E2  - A < L + l f L + l ) 

A<L+1 >L > - 0. 

RETURN ' . 

80  CONTINUE 
P - P/R 
Q - Q/R 

PREHULTIPLY. 
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00  90  J»LfN 
Z ■ A(Lf J) 

A(LfJ)  » P*Z  + Q*A<  L+l f J ) 

A < L+ 1 f J ) - P*A<L+1fJ>  - Q*Z 
90  CONTINUE 

POSTNULTIPLY. 

LI  - L+l 
DO  100  1*1 f LI 
Z =■  A(IfL) 

A ( I f L ) = P*Z  + Q*A<  I fL+1  ) 

A<I»L+1>  = P*A<IfL+1)  - Q*Z 
100  CONTINUE 

ACCUMULATE  THE  TRANSFORMATION  IN  V. 

DO  110  I*1fN 
Z - V(IfL) 

U< I fL ) = P*Z  + Q*V(IfL+1) 

V< I fL+1 ) » P*U<IfL+1)  - Q*Z 
110  CONTINUE 

A<L+1fL)  =*  0. 

El  =>  A ( L f L ) 

E2  = A<L+1fL+1) 

RETURN 

END 


SUBROUTINE  QRSTEP < Af VfP fQ fR f NL fNUf NfNA fNV ) 

INTEGER  NfNAfNLfNUfNV 
REAL  A<NAfN) fPfQfRfU<NUfN) 

QRSTEP  PERFORMS  ONE  IMPLICIT  OR  STEP  ON  THE 
UPPER  HESSENBERG  MATRIX  A.  THE  SHIFT  IS  DETERMINED 
BY  THE  NUMBERS  PfQf  AND  Rf  AND  THE  STEP  IS  APPLIED  TO 
ROUS  AND  COLUMNS  NL  THROUGH  NU.  THE  TRANSFORMATIONS 
ARE  ACCUMULATED  IN  V.  THE  PARAMETERS  IN  THE  CALLING 
SEQUENCE  ARE  (STARRED  APRAMETERS  ARE  ALTERED  BY  THE 
SUBROUTINE) 

*A  THE  UPPER  HESSENBERG  MATRIX  ON  UHICH  THE 

QR  STEP  IS  TO  BE  PERFORMED. 

*V  THE  ARRAY  IN  UHICH  THE  TRANSFORMATIONS 

ARE  TO  BE  ACCUMULATED 

*P  PARAMETERS  THAT  DETERMINE  THE  SHIFT. 

«Q 

*R 

NL  THE  LOUER  LIMIT  OF  THE  STEP. 

NU  THE  UPPER  LIMIT  OF  THE  STEP. 

N THE  ORDER  OF  THE  MATRIX  A. 

NA  THE  FIRST  DIMENSION  OF  THE  ARRAY  A. 

NV  THE  FIRST  DIMENSION  OF  THE  ARRAY  V. 

INTERNAL  VARIABLES. 

INTEGER  IfJfKfNL2fNL3fNUM1 
REAL  SfXfYfZ 
LOGICAL  LAST 

NL2  - NL+2 
DO  10  I-NL2fNU 
A< I f 1-2 ) » 0. 

10  CONTINUE 

IF(NL2  .EQ.  NU)  GO  TO  30 
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NL3  - NL+3 
00  20  I=NL3.NU 
A < I » l -3 ) - 0. 

20  CONTINUE 
30  CONTINUE 
NUH1  - NU-1 
DO  130  K»NL»NUN1 

DETERMINE  THE  TRANSFORMATION. 

LAST  = K .EQ.  NUM1 
IF(K  .EO.  NL)  GO  TO  40 
P = A<  K r K-l ) 

Q =*  A(K+1»K-1) 

R » 0. 

IF< .NOT. LAST)  R » A(K+2»K-1> 

X » ADS(P)  + ADS ( Q ) + ABS< R ) 

IF < X .EQ.  0.)  GO  TO  130 
P = P/X 
0 » Q/X 
R =*  R/X 
CONTINUE 
S » SORT  < P**2  + 

IF<P  .LT.  0.)  S 
IF< K .EQ.  NL)  GO  TO  50 
A<Kf K-l ) » -S*X 
GO  TO  60 
CONTINUE 

IF  (NL  .NE.  1)  A<  K rK-1 ) » -A(K»K-1> 
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40 


50 


60 


70 


80 


90 


100 


0**2  + R**2> 
* -S 


CONTINUE 
P = P + S 
P/S 
Q/S 
R/S 
Q/P 
R/P 


1 


PREMULTIPLT. 

DO  80  J»Kf N 

P = A < K f J ' + Q*A<K+1 • J ) 
IF(LAST)  GO  TO  70 
P = P + R*A<  K+2  f J) 

A< K+2 . J)  * A(K+2fJ>  - P*Z 
CONTINUE 

A<K+1 » J)  - P*Y 
A(KfJ)  - P*X 


A(K+1 » J) 
A(ICtJ) 
CONTINUE 


POSTMULTIPLY. 

J ■ MIN0(K+3»NU) 

DO  100  1*1 tj 

P - X*A< I f K ) + Y*A< I »K+1 ) 

IF (LAST)  GO  TO  90 
P - P + Z*A(I»K+2> 

A< I .K+2)  - A< I »K+2 ) - P*R 
CONTINUE 

A(  I »K  + 1 ) - A<I'K+1)  - P*Q 
A( I >K ) - A(I.K)  - P 
CONTINUE 

ACCUMULATE  THE  TRANSFORMATION  IN  U. 
DO  120  I-lfN 

P - /*V(IfK)  ♦ Y*V(I»K+1> 


. . 


SUBROUT I NE  ORTHES < NH t N » LOU , I GH , A . ORT ) 

INTEGER  IrJvHrNrIIr JJ r LA t HP t NM t IGH r KP1 r LOU 
REAL  A < NH  » N ) » ORT  < I GH  > 

REAL  F rGr Hr  SCALE 
LA  * IGH  - 1 
KP1  » LOW  + 1 
IF (LA  .LT.  KP1>  GO  TO  200 
DO  180  H-KPIfLA 
H = 0. 

ORT <H)  =*  0. 

SCALE  = 0. 

DO  90  1-1 r IGH  , 

90  SCALE  = SCALE  + ABS ( A< I r H-l ) ) 

IF ( SCALE  .EQ.  0.)  GO  TO  180 
HP  » H + IGH 
DO  100  II»HrIGH 
I - HP  - II 

ORT ( I ) - A < I f M- 1 ) /SCALE 
H =■  H + ORT  ( I ) *ORT  < I ) 

100  CONTINUE 

G » -SIGN (SORT (H) rORT <h) ) 

H » H-  ORT(H)*G 
ORT(M)  ■ ORT(H)  - G 
DO  130  J»HrN 
F - 0. 

DO  110  II»Hr IGH 
I * HP  - II 
F » F + 0RT<I)*A<IfJ> 

110  CONTINUE 
F - F/H 

DO  120  I-HfIGH 

A(IfJ)  » A< I r J)  - F*ORT(I) 

120  CONTINUE 

130  CONTINUE 

DO  160  I • 1 » I GH 
F ■ 0. 

DO  140  JJ-HrIGH 
J - HP  - JJ 
F - F + ORT(J)*A<I»J) 

140  CONTINUE 

F - F/H 

DO  ISO  J“H » I GH  • 

A<I»J)  - A( I » J)  - F*ORT(J> 

ISO  CONTINUE 

160  CONTINUE 

ORT (H)  - SCALEtORT ( H > 

A(H.H-l)  - SCALE6G 
180  CONTINUE 
200  RETURN 
END 


SUBROUTINE  ORTRAN<NM»NrLOMr IGHr ArORT fZ) 
INTEGER  I » JfNfKLfHHf HPrNHf IGHfLOWfHPI 
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REAL  A<NH » IGH ) fORT ( IGH ) r Z< NH  »N  > 
REAL  GfH 
00  80  I»1  fN 
DO  60  J-lfN 
Z(I»J)  ■ 0. 

60  CONTINUE 

Z(I»I)  « 1. 

80  CONTINUE 

KL  = IGH  - LOU  - 1 
IF <KL  .LT.  1)  GO  TO  200 
00  140  HH-lrKL 
HP  - IGH  - HM 
H ■ A(HPfHP-1)*0RT(HP> 

IF<H  .EQ.  0.)  GO  TO  140 
MP1  - HP+1 
DO  100  I»HP1 » IGH 
ORT(I)  » A< I fHP-1 ) 

100  CONTINUE 

60  130  J»MP  t IGH 
G = 0. 

DO  110  I=HPfIGH 

G - G + ORT  < I ) *Z ( I f J > 

110  CONTINUE 

G » G/H 

DO  120  I”MPf IGH 

Z<IfJ)  « Z(IfJ)  + G*ORT(I) 
120  CONTINUE 

130  CONTINUE 
140  CONTINUE 
200  RETURN 
END 


